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Abstract 



One of the major challenges in modern cosmology is to understand the reionization 
history of the Universe. This is directly related to galaxy formation and the formation 
of the first luminous objects. Observations of redshifted 21-cm radiation from neutral 
hydrogen (HI) is probably the most promising future probe of reionization. Several 
approaches have been proposed to extract information about the epoch of reionization 
from the data which is expected to come in near future. 

The most discussed approach has been to study the global statistical properties of 
the reionization HI 21-cm. We develop the formalism to calculate the Multi-frequency 
Angular Power Spectrum (MAPS) and quantify the statistics of the HI signal as a 
joint function of the angular multipole / and frequency separation Au. We adopt a 
simple model for the HI distribution which incorporates patchy reionization and use it 
to study the signatures of ionized bubbles on MAPS. We also study the implications of 
the foreground subtraction. 

This thesis also investigates the possibility of detecting ionized bubbles around in- 
dividual sources through radio interferometric observations of redshifted HI 21-cm ra- 
diation. We present a visibility based matched filter technique to optimally combine 
the signal from an ionized bubble and minimize the noise and foreground contributions. 
The formalism makes definite predictions on the ability to detect an ionized bubble or 
conclusively rule out its presence within a radio map. Results are presented for the 
GMRT and the MWA. Using simulated HI maps we analyzed the impact of HI fluctua- 
tions outside the bubble on its detectability. Various other issues such as (i) bubble size 
determination (ii) blind search for bubbles, (iii) optimum redshift for bubble detection 
are also discussed. 

Key words: cosmology: theory, cosmology: diffuse radiation, cosmology: large-scale 
structure of universe, Methods: data analysis 
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1 Introduction 



Understanding the evolutionary history of the Universe is one of the major goals in 
modern cosmology- Cosmic Microwave Background Radiation (CMBR) observations 
(COBBj, WM AP[J) give a picture of the early Universe (only ~ 370, 000 years after the 
Big Bang). During the first ~ 100000 years, the universe was a fully ionized plasma 
with a strong coupling caused by the Thomson scattering between photons and elec- 
trons. Because of adiabatic expansion, the temperature of the Universe dropped down 
to few thousand Kelvin (~ 3000 K) at redshit z ~ 1100 and the protons and electrons 
combined for the first time to form neutral Hydrogen (HI) atoms. The scattering of 
photons reduced and they decoupled from baryonic matter. After this the photons were 
mostly undisturbed except that the expansion of the Universe redshifted them into the 
microwave at present. This relic background radiation is known as the CMBR. After 
the recombination, the Universe became almost neutral with 75% and 25% in weight of 
total baryonic matter was in the form of HI and neutral Helium respectively (neutral 
Helium formed earlier than the HI). On the other hand observations of Ly alpha forest 
in quasar (QSO) absorption spectra sho w the diffu s e Hyd rogen gas in the Universe to 



be completely ionized at redshifts z < 5 (IFan et al 



(|2002h ) 



*http://lambda. gsfc.nasa.gov/product/cobe/ 
thttp://map. gsfc.nasa.gov/ 
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1 Introduction 



1.1 The Epoch of Reionization 



The above two observations suggest that the HI of the Universe was ionized sometime in 
the redshift range z = 1000 to 5. The period when the HI was ionized is called the Epoch 
of Reionization (EoR). The EoR is one of the least known chapters in the evolutionary 
history of the Universe. Its exact timing, duration, nature of the reionizing sources, 
their relative contribution to the reionization, large scale distribution and evolution of 
HI are highly unknown. Again, once the first reionizing sources were formed, their 
various feedback mechanisms such as mass deposition, energy injection and emitted 
radiation deeply affect subsequent galaxy formation and influence the evolution of the 
Intergalactic Medium (IGM). The epoch of reionization, therefore, can be considered 



as a comp 
review see 



icated era which invo 



v es a l arge number of interconnected processes (for a 



Choudhury &: Ferraral (120061 ) ). The study rei o nizati o n has been at t he fore 



font of research over the last few years (jBarkana fc Loebl (120011 ); 



Barkanal (120061 )) 



Though the EoR is yet to be observed in detail, theorists have proposed possible 
pictures of how the reionization took place. The weak density perturbations which 
were generated during inflation era grew through gravitational instability and lead to 
overdense regions. The first generation of galaxies formed at redshift z ~ 20 in these 
overdense regions. The gas in these galaxies cooled by molecular cooling and fragmented. 
Then the first generation of stars which are believed to be massive (~ 100M o ) and metal 
free were created. Enormous amount of radiation produced by these stars ionized the 
surrounding IGM. Then new generation of galaxies and stars formed. Ionized bubbles 
thus grow and filled the entire space. In another scenario, black hole were created at the 
centers of galaxies. Enormous amount of x-ray radiation ionized the IGM. 

Currently two types of experiments give information about the EoR. First, the study 
of Lyman-a line absorption in the high redshift QSO spectra has been be used to probe 



2 



1.1 The Epoch of Reionization 



the ionization sta te and the HI distributi on at high redshifts. The analysis of Gunn- 



Peterson troughs (jGunn fc Peterson! (119651 )) in the high re dshift QSO absorpt i on spectra 



suggests that reionizatio n finished around redshift z ~ 6 ( IBecker et al. 



2001 



Fan et al. 



2002 : 



White et al 



20031 ). The study of the dark gap distribution and its evolution in 
QSO spectra has been shown as an efficient and independent probe of the reionization. 
his analysis puts an upper limit on the HI fraction xhi < 0.36 at redshift z = 6.3 



( Gallerani et al. 



20081 ) . Measurements of the siz es of HII r e gions around high redshift 



(120061 ) ) . Another indepen- 



QSOs are also consistent with the above results (IFan et al. 
dent constraint comes from CMBR observations. The CMBR photons scatter off free 
electrons (produced during reionization) which results in the suppression of the intrinsic 
temperature and polarization anisotropies on angular scale below the horizon at EoR. 
At the same time a polarization signal is generated at large angular scale. The ampli- 
tude and the position of the peak in the polarization angular power spectrum depend 
on the reionization redshift. From recent WMAP measurements of the electron scat- 



tering optical de pth r P from the temperature-po 



power spectrum (IPage et al. 



2007 



Dunkley et al. 



arization and polarization-p olarization 



2008 



Komatsu et al 



2008|) imply that 



reionization started before z ~ 10. It thus seems from these two experimental results 



that the EoR is an extended and complex process which occurred over a redshift 



6—15 (jChoudhury fc Ferrara 



2006a : 



Fan. Carilli &: Keating 



2006 



Alvarez et al. 



rang e 



20061 ) 



However, there exist limitations of using these observations to study the details of reion- 
ization. The Lyman-a absorption feature (similarly for Lyman-/? and Lyman-7) in QSO 
spectra is not sensitive to the higher neutral fraction of hydrogen (> 0.1). Therefore, 
these observations can not be used to probe the reionization at its earlier stages. The 
CMBR experiments are sensitive only to the integrated history of the EoR and it may not 
be useful to study the progress of reionization with redshift. In fact, it has been shown 
that the CMBR polarization power spectrum is weakly dependent on the details of the 
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reionization his tory (jKaplinghat et al 



2003 



Hu fe Holder 



2003 



Haiman Holder 



2003 



Colombo 



20041 ). though weak constraints could be obtained from upcoming experiments 
such as PLANCKE. 

Several other probes have been discussed in the literature to unveil the reionization 
history of the Universe. According to theory we expect QSOs and massive galaxies to 
form in highly overdense regions. Moreover first generation stars which are expected 
to be massive and short lived would produce gamma ray bursts. Observations of high 
redshift QSOs, galaxies and gamma ray bursts could in principle provide a great deal of 
information. Though currently available experimental sensitivity is not sufficient enough 
to detect those, future space based experiments like JWSTlj will have enough sensitivity 
to detect those objects. 



1.2 21-cm Tomography 

The interaction between the spins of the proton and the electron in a hydrogen atom in 
its ground state gives rise to two hyperfine states, i.e., the triplet states of parallel spins 
and the singlet state of anti-parallel spins. The triplet state has higher energy than the 
singlet state. When a hydrogen atom jumps from the triplet to the singlet state it emits 
a photon with the wavelength of 21-cm. 

Observations of the redshifted 21 cm line from the EoR i s perhaps one of the mo st 



Furlanetto et al. 



feood)). 



promising tools for studying the EoR (for recent review see 
The advantage of this probe lies in the fact that the EoR can be probed at any desired 
redshift by appropriately tuning the observation frequency. Since HI is distributed all 
over space these observations have the potential to probe the large scale distribution of 



^http: / / www.rssd.esa.int /Planck / 
§http: / / www.jwst.nasa.gov/ 
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HI. A wealth of information about the EoR can be extracted from these observations. 
Unlike the observations of QSO absorption spectra the redshifted 21 cm radiation does 
not suffer from saturation because the optical depth for 21 cm radiation is much less 
than unity. 

The possibility of o bserving 21 cm emission from t he cosmological stru cture formation 



was fir s t recognized by 



(119791 ) 



Scott fc Reesl (119901 ) and 



Sunvaev and Zeldovich ( 



1972) and 



ater st udied by 



Madau. Meiksin k, Reesl f . 19971 ) considering both emis- 



Hogan &: Rees 



sion and absorption against the CMBR. More recently, the effect of heating of the HI 
gas and its reionizatio n on 21 cm signal has d r awn great deal of attention and has 



been studied in detai 



Iliev et al 



2002 



2003 



Gnedin fc Ostri 



kcr 



Ciardi fc Madau 



1997 



2003 



Shaver et al 



1999 



Tozzi et al. 



Furlanetto. Sokasian &: Hernquist 



2000; 



2004 



2005; 



Mcquinn et al. 



Miralda-Escude!l2003:IChen fc Miralda-Escudelbooi ICoorav & Furlanetto 120051: ICoorav 



2006 



Seth 



2005 



Salvaterra et al 



2005 



Carilli 



20061 ). 



There could be several approaches in interpreting the data which is expected to 
come in coming years. Measurements of 21-cm signal in emission averaged over large 
area of sky would p r ovide a direct probe of t he ey olution of neutral fraction with 



redshift (IShaver et al. 



(1999) 



Gnedin fc Shaverl (120041 ) ) . Telescopes are being set up 



to measure this average signal (e.g, Compact Reionization Experiments (CORE) at 
Australia Telescope National Facility, Experiments to Detect Global EOR Signature 
(EDGES) at MIT Haystack Observatory). Measurements of the HI power spectrum and 
individual ionized bubble detection are other two major approaches which are discussed 
in the following subsections. 
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1.2.1 Statistical analysis of the 21-cm signal 



The most discussed approach has been to study the global statistical properties of the 
HI distribution through quantities like the power spectrum. The precise measurements 
of the HI 21-cm fluctuations in terms of their multifrequency angular power spectrum 
would provide a wealth of information of the size, spatial distribution and evolution 
of the ionized regions. This would also help us understand the effect of reionization 
on the structure formation, radiative feedback mechanisms in star- forming zones, the 



2001 



Bharadwaj fc Pandey 



2003 



(Bharadwai. Nath & Sethi 


2001; 


Bharadwai &; Sethi 


Bharadwaj & Srikant 


2004 


). A similar formalism can 



also be applied at high redshifts to probe reionization and also the pre-reionization era 



( Zaldarriaea. Furlanetto & t 


ernauist 


2004; 


Furlanetto. Zaldarriaea & Hernauist 


2004b: 


Bharadwa 


& Ali 


2004 


2005 




B 


raradwai & Pandev 


2005: 


Ali. Bharadwai and Pandev 


2005. 


2006; 


Loeb & Zaldarriaga 


2004; 


He et al. 


2004 


). It is expected that the radiation 



from first generation luminous objects changes the character of the 21 cm sky completely. 
During this epoch, an unique signature of ionized regions will be imprinted on the 
redshifted 21 cm signal that manifests the processes for the ionizing radiations and that 
evolves with redshift as reionization proceeds. Chapter 2 presents the multifrequency 
angular power spectrum (MAPS) of the epoch of reionization 21-cm signal as a joint 
function of the angular multipole I and frequency separation Av. This studies the 
signature of ionized regions (bubbles) on the MAPS and its implication for separating 
foregrounds from the signal. 
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1.2.2 Detecting individual ionized bubbles in 21-cm maps 



It is believed that the ionizing radiation from QSOs and the stars within galaxies 
reionize the surrounding neutral IGM. The initial fram ework for the growth of HII 



regions around ind i vidua l galaxies have been developed by lArons fc Wingert 



(11972) and 



Shapiro fc Giroux! (119871 ). Later various types of sophisticated models for the growth 



of HII regions are p rescribed and used to sharpen our understanding about inhomo 



geneous reionizaton (IFurlanetto. Zaldarriaga fc Hernquist 



2004a 



Kramer et al. 



Cohn fc Chang 



2007 



20061 ). A different, complimentary approach would be to directly ob- 
serve the individual ionized regions around luminous sources (stars/QSOs). The issue of 
detecting these bubbles in radio-interfero metric observations of redshifted HI 21 cm radi- 
ation has been drawing considerable attention. The detection of individual ionized bub- 
bles would be a direct probe of the reionization process. It has also been proposed that 
such observations will probe the properties of the ionizing sources and the evolution of the 



surrounding IGM f 1 


vVvithe & Loeb. 


2004a: 


Wvithe. Loeb &; Barnes 


200 


5 


Kohler et al.. 


2005; 


Mase 


li et al. 


2007: 


Alvarez & Abel 


2007: 


Geil & Wvithe 




2008 




Wvithe 




2008: 


Geil et al. 


> 


2008 


). Observations of individual ionized bubbles would complement the 



study of reionization through the power spectrum of HI brightness temperature fluctu- 
ations. 

Nearly all the above mentioned work on detecting ionized regions consider the con- 
trast between the ionized regions and the neutral IGM in images of redshifted HI 21 cm 
radiation. The HI signal is expected to be only a small contribution buried deep in the 
emission from other astrophysical sources (foregrounds) and in the system noise. Chap- 
ter 3 introduces a matched filter to optimally combine the entire signal of an ionized 
bubble while minimizing the noise and foreground contributions. This technique uses 
the visibilities which are the fundamental quantity measured in radio-interferometric ob- 
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servations. Using visibilities has an advantage over the image based techniques because 
the system noise contribution in different visibilities is independent whereas the noise 
in different pixels of a radio- interferometric images is not. Chapter 4 presents simula- 
tion results for bubble detection in 21-cm maps and also studies the size and position 
determination of ionized bubbles. 



1.3 The Radio- Interferometric Experiments 

On the experimental side several low frequency radio experiments are either functional 
or being set up. This motivated us to study the expected redshifted 21 cm background 
from the EoR and possibility to detect it. Two different observational strategies will 
be followed. The first approach is to measure the global evolution of mean HI signal 
with redshift, and second is to measure large scale distribution of HI through the power 
spectrum measurements and detect the HII regions. In principle both approaches would 
provide a wealth of information about the reionization. The global signature experi- 
ments use a single, small antenna which provides a large field of view (FoV). These will 
measure the HI 21-cm signal in emission averaged over a large area of the sky. Two 
main experiments, namely the Cosmological Reionization Experiments (CORE) at the 
Australian Telescope National Facility, and the Experiment to Detect the Global EoR 
signature (EDGES), at the MIT Haystack Observatory Q are underway in this direction. 

The majority of the recent or upcoming radio-interferometric experiments are aimed 
at measuring the HI 21-cm signal statistically. Individual ionized wholes (bubbles) 
can also be detected using interferometric obs ervations of HI. The Giant Metre- Wave 



Radio Telescope (GMRTlI 



Swarup et al. 



1991) is already functioning at several bands 



"http: / / www. hay s t ack .mit.edu 
II http: / / www.gmrt.ncra.tifr.res.in 
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in the frequency range 150-1420 MHz and can potentially detect the 21 cm signal at 
high redshifts. In addition, construction of other low-frequency experiments such as the 
Murchison Widefield Array (MWA), LOw Frequency ARray (LOFAFo), 21 Centemeter 
Array (21CMA[2I), Precision Array to Probe Epoch of Reionization (PAPER), Square 
Kilometer Array(SKAEj) has raised the possibility to detect 21 cm signal from very high 
redshifts. The above first generation experiments will probably start their operation at 
the end of the decade. 



1.4 Challenges to Overcome 



Although the redshifted 21-cm line can provide enormous amount of information, its 
detection is going to be a huge challenge. The signal is expected to be highly con- 
taminated by foreground radio emission. Potential sources for these foregrounds include 
synchrotron and free-free emission from our Galaxy and external gala xies, low-frequency 



radio point sources and free 



DiMatteo. et al. 



2002 . 



2004 



ree emission from electrons in t 



Oh k Mack 



2003 



Gleser et al. 



2008 



re IGM (Shaver et al. 



Ali et al 



1999 



20081 ). Con- 



tributions from astrophysica 



foregrounds are expected to b e several order of magnitude 



stronger than the HI signal (jSantos. Cooray fc Knox! (120051 )) 



However, there have been various proposals for tackling the foregrounds, the most 
promising being the application of multi-frequency observations. The foregrounds are 
expected to have a continuum spectra, and the contribution at two different frequencies 
separated by Av ~ 1 MHz are expected to be highly correlated. The HI signal, on 
the other hand, is expected to be uncorrelated at such a frequency separation and this 
holds the promise of allowing us to separate the signal from the foregrounds. It has been 



**http://www. lofar.org/ 

tthttp://web.phys. cmu.edu/~past/ 
"http://www.skatelescope.org/ 
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proposed that multi- 
the fo reground (e.g. 



frequency analysis of the radio signal can be useful in separating out 



Zaldarriaga. Furlanetto fc Hernquist 



2004 



Santos. Cooray fc Knox 



20051 ). A n alternate approa ch is to subtract a best fit continuum spectra along each line 



of sight (jWang et al 



20061 ) and then determine the power spectrum. This is expected 



to be an effective fo reground subtraction method in data with very low noise levels. 



Morales et al. 



(120061 ) have discussed the complementarity of different foreground removal 
techniques and the implications for array design and the analysis of reionization data. 

The system noise in all low frequency radio experiments relevant for the reionization 
is dominated by the sky contribution T sky . We also expect T sky to vary depending on 
whether the source is in the galactic plane or away from it. This is expected to be an 
independent Gaussian random variable and can be reduced by increasing observation 
time. The observation time required to detect the HI signal is an important issue. 

Man-made radio frequency interference (RFI) is a growing problem in all earth-based 
radio astronomy. Signal from television, FM radio, satellites, mobile communication, 
electric spark etc. all fall in the same frequency band as the redshifted 21-cm reion- 
ization signal from the reionization. These are expected to be much stronger than the 
expected 21cm signal, and it is necessary to quantify and characterize the RFI. It has 
been suggested to construct EoR experim ents at remote sites which are expected to have 



low RFI. Recently 



Bowman et al 



(120071 ) have characterized the RFI for the MWA site 
on the frequency range 80 to 300MHz. They find an excellent RFI environment except 
for a few channels which are dominated by satellite communication signal. There are 
several methods and techniq ues being developed to identify, characte rize, and ultimately 



subtract interfering signals (IFridman fc Baan 



2001 



Ellineson 



2005). Telescope design 



2000). The effect of 



also plays a role in mitigating interfering signals (ILeshem et al. 
polarization leakage is another issue which needs to be investigated in detail. This could 
cause polarization structures on the sky to appear as frequency dependent ripples in 
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the foregrounds intensity . This could be particularly severe for the MWA which has 
a wide field of view. Radio recombinations lines could be a significant contaminant in 
the low frequency radio observations. Unfortunately due to lack of observations we have 



l ittle knowledge o 



(Morales &: Hewitt 



the i mpact of these lines in the epoch of reionization observation 



2004]). Refraction index of Earth's atmosphere varies significantly 
in space as well as in time in the frequency band relevant for EoR experiments. This 
creates significant calibration and imaging proble ms that must be solved in order to re - 



1986|). 



liably clean the strong foreground contamination ([Thompson. Moran fc Swensonl . 
This is another important issue related to system noise. The FoV of the individual an- 
tenna and the baseline distribution change with observing wavelength and if these are 



neglected they could cause severe problem in extracting the signal. 



1.5 Outline of the Thesis 



We give an outline of the rest of the thesis. 



Chapter 2 calculates the Multi-frequency Angular Power Spectrum (MAPS) to 
quantify the statistics of the HI signal as a joint function of the angular multipole I and 
frequency separation Au. Assuming a small portion of a spherical sk y as a flat-sky we 



develo p formulae for MAPS, including the effect of peculiar velocities ( jBharadwaj fc Ali 



(120041 ) ). The flat sky approximation is found to be a good representation over the 
angular scales of interest. The final expression is very simple to calculate and interpret 
in comparison to the formulae obtained using full spherical sky. We adopt a simple model 
for the HI distribution which incorporates patchy reionization and use it to study the 
signatures of ionized bubbles on MAPS. We also study the implications of the foreground 
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subtraction. 

Chapter 3 investigates the possibility of detecting individual ionized regions (bub- 
bles) in radio-interferometric observations of HI 21 cm radiation. We develop a visibility 
based formalism that uses a matched filter to optimally combine the entire signal from 
a bubble while minimizing the noise and foreground contributions. The method makes 
definite predictions on the ability to detect an ionized bubble or conclusively rule out 
its presence within a radio map. We make predictions for the GMRT and the MWA at 
a frequency of 150 MHz (corresponding to a redshift of 8.5). 

Chapter 4 studies the impact of the HI fluctuations outside the bubble that we 
are trying to detect on the detectibility of the bubble in 21-cm maps. We use simulated 
HI maps which incorporates the patchy reionization scenarios and investigate the re- 
strictions imposed by the HI fluctuations on bubble detection. We validate the matched 
filter technique presented in Chapter 3 through simulation of bubble detection. We also 
use the simulations to determine the accuracy to which the GMRT and the MWA will 
be able to determine the size and position of an ionized bubble, and test if this is limited 
due to the presence of HI fluctuations. 

In Chapter 5 we estimate the optimum redshift for detecting ionized bubbles in 21- 
cm maps for different reionization scenarios. We investigate the situations under which 
bubbles can be detected. Results are also presented in terms of the scaling relations. 
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2 Multi-frequency Angular Power 



Spectrum of 21 cm Signal 



2.1 Introduction 

In this Chapter, we develop the formalism to calculate the multi-frequency angular 
power spectrum (hereafter MAPS) which can be used to analyse the 21 cm signal from 
HI both in emission and absorption against the CMB. We restrict our attention to HI 
emission which is the situation of interest for the epoch of reionization. In our formalism, 
we consider the effect of red shift-space distort i ons w hich has been ignored in many of 



earlier works. As noted by 



Bharadwaj &: Alii (120041 ). this is an important effect and 



can enhance the mean signal by 50% or more and the effect is expected to be most 
pronounced in the multi-frequency analysis. We next use the flat sky approximation 
to develop a much simpler expression of MAPS which is much easier to calculate and 
interpret than the angular power spectrum written in terms of the spherical Bessel 



This chapter is adapted from the paper " The multi-frequency angular power s p ectrum 
of the epoch of reionization 21 cm signal" bv \Datta. Choudhurv & Bharadwa-i \20Q'A ). 
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functions. We adopt a simple model for the HI distribution (IBharadwaj fc Alii . 120051 ) 
which incorporates patchy reionization and use it to predict the expected signal and 
study its multi-frequency properties. The model allows us to vary properties like the 
size of the ionized regions and their bias relative to the dark matter. We use MAPS to 
analyze the imprint of these features on the HI signal and discuss their implication for 
future HI observations. 

As noted earlier, the HI signal at two different frequencie s separated by Ay is ex - 



Bharadwai fc Alii (120051 ) 



pected to become uncorrelated as Au is increases. As noted in 
the value of Au beyond which the signal ceases to be correlated depends on the an- 
gular scales being observed and it is < 1 MHz in most situations of interest. A prior 
estimate of the multi-frequency behavior is extremely important when planning HI ob- 
servations. The width of the individual frequency channels sets the frequency resolution 
over which the signal is averaged. This should be chosen sufficiently small so that the 
signal remains correlated over the channel width. Choosing a frequency channel which 
is too wide would end up averaging uncorrelated HI signal which would wash out various 
important features in the signal, and also lead to a de gradation in the signal to noise 



2005J) 



ratio. In this context we also note that an earlier work (jSantos. Cooray fc Knox 
assumed individual frequency channels 1 MHz wide and smoothed the signal with this 
before performing the multi-frequency analysis. This, as we have already noted and shall 
study in detail in this Chapter, is considerably larger than the Au where the signal is 
uncorrelated and hence is not the optimal strategy for the analysis. We avoid such a pit- 
fall by not incorporating the finite frequency resolution of any realistic HI observations. 
It is assumed that the analysis be used to determine the optimal frequency channel 
width for future HI observations. Further, it is quite straightforward to introduce a 
finite frequency window into our result through a convolution. 

The outline of this Chapter is as follows. In Section 2.2 we present the theoretical 
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formalism for calculating MAPS of the expected 21 cm signal considering the effect of 
HI peculiar velocity. The calculation in the full-sky and the flat-sky approximation are 
both presented with the details being given in separate Appendices. Section 2.2.1 defines 
various components of the HI power spectrum and Section 2.2.4 presents to models for 
the HI distribution. We use these models when making predictions for the expected HI 
signal. We present our results in Section 2.3 and also summarize our findings. In Section 
2.4 we discuss the implications for extracting the signal from the foregrounds. 



2.2 Theoretical Formalism 



2.2.1 The HI power spectrum 



The aim of this Section is to set up the notation and calculate the angular correla- 



tion Ci for the 2 



Bharadwaj fc Ali 



cm b rightness temperature fluctuations. It is now well known (e.g.. 



20051 ) that the excess brightness temperature observed at a frequency 



v along a direction n is given by 



T(u, n) = T(z) rjm (z, hr u ) 



(2.1) 



where the frequency of observation is related to the redshift by v = 1420/(1 + z) MHz. 
We consider a flat Universe (k = 0) in which the comoving distance r u can be written as 



(2.2) 



The mean background excess brightness temperature T(z) at redshift z is written as 



T(z) » 25mK 



0.15 ffl b h 2 \ fl-Y\ 1 + z 



Vl m h 2 V 0.022/ V 0.76 



10 



(2.3) 
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where Y « 0.24 is the helium mass fraction and all other symbols have usual meaning. In 

1 /2 

the above relation, it has been assumed that the Hubble parameter H\z) ~ HoQm (1 + 
z) 3 / 2 , which is a good approximation for most cosmo logical models at z > 3. The quan- 



tity r ?HT is known as the "21 cm radiation efficiency in redshift space" (IBharadwaj fc Ali 



20051 ) and can be written in terms of the mean neutral hydrogen fraction xhi and the 



fluctuation in neutral hydrogen density field <5hi as 



77Hi(2,nr„ 



x m (z)[l + 5 m (z,hr u )] (^-Tf) 
(1 + z) dv(z, nr„ 



H(z) 



(2.4) 



where T 7 and T s are the temperature of the CMB and the spin temperature of the gas 
respectively. The term in the square bracket arises from the coherent components of 
the HI peculiar velocities. In the above derivation it is assumed that the term (1 + 
z)dv(z, nr v )/dr v is small compared with H(z) which is a reasonable assumption for the 
scale of our interest. 

At this stage, it is useful to make a set of assumptions which will simplify our anal- 
ysis: (i) We assume that T s ^> T 7 , which corresponds to the scenario where the spin 
temperature T s and the gas kinetic tempera ture are strongly coupled eithe r through 



strong Lya scattering or collisional coupling fMa dau. Meiksin fc Rees 



the couplings are expected to be patchy (IHiggins fc Meiksin 



. 1997ft . Though 



20091 ) the assumption is 



reasonable throughout the IGM soon after the formation of first sources of radiation, 
(ii) We assume that the HI peculiar velocity field is determined by the dark matter 
fluctuations, which is reasonable as the peculiar velocities mostly trace the dark matter 
potential wells. This assumption is valid for scales larger than the Jeans length scale 
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which are the scales of our interest. We then have 

/A3L 

_^_ e -i fc ,,(k.n)~ Hi( ^ k) (25) 

where for k 7^ 

fj m (z, k) = x m (z) [A H1 (z, k) + (k • h) 2 A(z, k)] , (2.6) 

and Abi(z, k) and A(z, k) are the Fourier transform of the fluctuations in the HI and the 
dark matter densities respectively. Note that /(f2 m ), which relates peculiar velocities to 
the dark matter, has been assumed to have a value f(tt m ) = 1 which is reasonable at 
the high z of interest here. 

For future use, we define the relevant three dimensional (3D) power spectra 

(A(z,k)A*(z,k')> = (2n) 3 5 D (k-k')P(z,k) 
(A m (z,k)A* m (z,k')} = (27rf5 D (k-k')P A ,Jz,k) 
(A(z,k)A* m (z,k')) = (27rf5 D (k-k')P Am (z,k) (2.7) 



where P(z, k) and P A ^(z, k) are the power spectra of the fluctuations in the dark matter 
and the HI densities respectively, while Pa h1 (z, k) is the cross-correlation between the 
two. 
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2.2.2 The multi-frequency angular power spectrum (MAPS) 



The multi-frequency angular power spectrum of 21 cm brightness temperature fluctua- 
tions at two different frequencies v\ and z/ 2 is defined as 

G\{y\, v 2 ) = (ai m (ui) a* m (u 2 )} . (2.8) 

In our entire analysis V\ and v 2 are assumed to differ by only a small amount Av <C v\, 
and it is convenient to introduce the notation 

Ci{Au) = Ci{u,u + Au) (2.9) 

where we do not explicitly show the frequency v whose value will be clear from the 
context. Further, wherever possible, we shall not explicitly show the z dependence of 
various quantities like T, xm, P{k) etc., and it is to be understood that these are to be 
evaluated at the appropriate redshift determined by v. 

The spherical harmonic moment of T(u, n) are defined as 

ai m {v) = J 'dnyj^(n) T[v,n) 

= TjdnY^h) J i^Hi (k) e- ifc ^ fi ) . 

(2.10) 

Putting the expression (12.61) for fjm (k) in the above equation, one can explicitly calculate 
the MAPS in terms of the three dimensional power spectra defined earlier. We give the 
details of the calculation in Appendix |A] and present only the final expression for the 
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angular power spectrum at a frequency v 

2f 2 xl 



<yri r°° r 
Q{Au) = S / k 2 dk Ukr v ) 3l (kr U2 )P A2 (k) 

7T Jo L HI 

-{jl(kr v )ji(kr U2 ) + ]i{kr U2 )j"(kr u )}P Am (k) 
+fi'(kr u )j' l '(kr V2 )P(k)] 



(2.11) 



Here j"{x) = -^ji(x) and we have used the notation r„ 2 = r v + r' v Au with 

r^ s g^ = -£( 1+ /) a . (2.12) 
" dv u H{z) K } 

Note that equation (12. lip predicts Ci(Au) from the cosmological 21 cm HI signal to be 
real. 

With increasing Au, we expect the two spherical Bessel functions ji(kr Ul ) and ji(kr U2 ) 
to oscillate out of phase. As a consequence the value of Ci(Au) is expected to fall in- 
creasing Au. We quantify this through a dimensionless frequency decorrelation function 
defined as the ratio 

For a fixed multipole I, this fall in this function with increasing Ai> essentially measures 
how quickly features at the angular scale 9 ~ tt/1 in the 21 cm HI maps at two different 
frequencies become uncorrelated. Note that < |k/(Az/)| < 1. 



2.2.3 Flat-sky approximation 

Radio interferometers have a finite field of view which is determined by the parameters 
of the individual elements in the array. For example, at 150 MHz this is around 3° for the 
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GMRT. In most cases of interest it suffices to consider only small angular scales which 
correspond to I 3> 1. For the currently favored set of flat ACDM models, a comoving 
length scale R at redshift z > 5 would roughly correspond to a multipole 

I « 3 x 10 4 (l - -±L=) ( 1 (2-14) 

Thus, for length scales of R < 100/i -1 Mpc at z ~ 10, one would be interested in 
multipoles I > 200. For such high values of I one can work in the flat-sky approximation. 

A small portion of the sky can be well approximated by a plane. The unit vector n 
towards the direction of observation can be decomposed as 

n = m + #; m ■ 9 = 0; | 9 |< 1 (2.15) 

where m is a vector towards the center of the field of view and 9 is a two-dimensional 
vector in the plane of the sky. It is then natural to define the two-dimensional Fourier 
transform of T(u, n) in the flat-sky as 

f(u,V) = [ d 2 9 e- 27TiU - e 2>, n) (2.16) 



where U, which corresponds to an inverse angular scale, is the Fourier space counter- 
part of 0. Using equations (12.11) and (12.51) . and the fact that for the flat-sky we can 
approximate k • n m k • m = k» we have 

f{u, U) = ^ J dk\\ e~ ik ^ fki(k\\m + 2nU/r u ) . (2.17) 
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It is useful to introduce the ff^i power spectrum P m defined as 

(%i(k)%i(k')> = (2vr) 3 4(k - k')Pffl(k) • (2.18) 
This is related to the other three power spectra introduced earlier through 

P ffl (k) =x 2 m (z)[P A2 Jk) + 2^P Am (k)+^P(k)} (2.19) 

I II 1 

where \i = m ■ k = k\\/k ( iBarkana fc Loebl 120051 ) . Note that the anisotropy of Pm(k) 

i.e., its /i-dependence arises from the peculiar velocities. 

The quantities calculated in the flat-sky approximation can be expressed in terms 
of their all-sky counterparts. The correspondence between the all-sky angular power 
spectra and its flat-sky approximation is given by 

(f{u 1 ,V)f*{u 2 ,TJ f )) = C^ u {u 1 ,u 2 ) 4 2) (U-U') (2.20) 

where 5p (U — U') is the two-dimensional Dirac-delta function. The details of the above 
calculation are presented in Appendix El Thus allows us to estimate the angular power 
spectrum C\ under the flat-sky approximation which has a much simpler expression 

ry2 POO 

Cf at (A^ = — / dfc|| cos(fc||<Az/) P ffl (k) (2.21) 

where the vector k has magnitude k = + P/r% ie. k has components fcy and l/r v 
along the line of sight and in the plane of the sky respectively. It is clear that the 
angular power spectrum Ci(Au) is calculated by summing over all Fourier modes k 
whose projection in the plane of the sky has a magnitude l/r v . We also see that Ci is 
determined by the power spectra only for modes k > l/r u . 
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The flat-sky angular power spectrum Cf at (0) is essentially the 2D power spectrum of 
the HI distribution on a plane at the distance r v from the observer, and for Av = equa- 
tion^ is just the relation between the 2D power spectrum and its 3D counterpart 



(jPeacock 



(11999I )). For Av ^ it is the cross-correlation of the 2D Fourier components 
of the HI distribution on two different planes, one at r v and another at r u+/ \ u . Any 
2D Fourier mode is calculated from its full 3D counterparts by projecting the 3D modes 
onto the plane where the 2D Fourier mode is being evaluated. The same set of 3D modes 
contribute with different phases when they are projected onto two different planes. This 
gives rise to cos(k\\r' u Av) in equation ( 12.21ft when the same 2D mode on two different 
planes are cross-correlated and this in turn causes the decorrelation of Cj (Av) with 
increasing Av. 

Testing the range of I over which the flat-sky approximation is valid, we find that 
for the typical HI power spectra Cf at (Az/) is in agreement with the full-sky Ci(Av) 
calculated using equation (12.111) at a level better than 1 per cent for angular modes 
I > 10. Since the integral in equation (12.211) is much simpler to compute, and more 
straightforward to interpret, we use the flat-sky approximation of C/ for our calculations 
in the rest of this Chapter. 

Note that equation (12.211) is very similar to the expression for the visibility correla- 



tions [equation (16) of 



Bharadwaj fc Ali 



20051 ] expected in radio interferometric observa- 



tions of redshifted 21 cm HI emission. The two relations differ only in a proportionality 
factor which incorporates the parameters of the telescope being used for the observation. 
This reflects the close relation between the visibility correlations, which are the directly 
measurable quantities in radio interferometry, and the C/s considered here. 
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2.2.4 Modeling the HI distribution 



The crucial quantities in calculating the angular correlation function are the three di- 
mensional power spectra P(k), PA m {k) an d Pa 2 {k). The form of the dark matter power 



spectrum P(k) is relatively well-established, pa rticularly within 



;he lin ear theory. We 



. 1997h ) 



shall be using the standard expression given by (IBunn fc White! 

The power spectrum of HI density fluctuations P^p. (k) and its cross-correlation with 
the dark matter fluctuations PA m (k) are both largely unknown, and determining these is 
one of the most important aims of the future redshifted 21 cm observations. A possible 
approach could be to implement some specific model for reionization, then attempt to 
predict the expected patchiness in the HI distribution and calculate the power spectra. 
Such an exercise is somewhat beyond the scope of this paper. The objective here is to 
quantify the angular power spectrum in terms of the physical attributes characteristic 
of the HI distribution at the epoch of reionization. To this end we adopt two simple 
models with a few parameters which capture the salient features of the HI distribution. 

The first model, which we shall denote as DM, assumes homogeneous reionization 
where the HI traces the dark matter , i.e., Am = A. This model does not introduce 
a characteristic length-scale in the HI distribution, and hence it serves as the fiducial 
model against which we can compare the predictions for patchy reionization. Under the 
standard scenario of reionization by UV sources, this is a valid assumption in very early 
stages of reionization when most of the IGM is neutral. However, this assumption could 
have another range of validity. This has to do with the scenarios where the dominant 
source of reionization are the exotic decaying particles, like neutrinos. In such case, 
there would be no bubbles associated with individual galaxies, rather the reionization 
proceeds in a homogeneous manner. In this model we have 



23 



2 Multi-lrequency Angular Power bpectrum of 21 cm bignal 



P A *Jk) = P Am (k) = P(k) (2.22) 
which we use in equation (12.111) to calculate Ci(Au). Alternately, we have 

Pw(k) =x 2 m (l + fi 2 ) 2 P(k) (2.23) 

which we can use in equation (I2.21j) to calculate C;(Az/) in the flat-sky approximation. 
This model has only one free parameter namely the mean neutral fraction xm- 

The second model, denoted as PR, incorporates patchy reionization. It is assumed 
that reionization occurs through the growth of completely ionized regions (bubbles) in 
the hydrogen distribution. The bubbles are assumed to be spheres, all with the same 
comoving radius R, their centers tracing the dark matter distribution with a possible 
bias b c . While in reality there will be a spread in the shapes and sizes of the ionized 
patches, we can consider R as being the characteristic size at any particular epoch. 
The distribution of the centers of the ionized regions basically incorporates the fact 
that the ionizing sources are expected to reside at the peaks of the dark matter density 
distribution and these are expected to be strongly clustered. For non-overlapping spheres 
the fraction of ionized volume is given by 

ArrR 3 

xnn = 1 - x m = —^—n HU (2.24) 

where nun is the mean comoving number density of ionized spheres and we use TZ to 
denote the ratio TZ = xhii/^hi- This model has been discussed in detail in Bharadwaj 
and Ali (2005), and we have 

A ffl (k) = [1 - b c KW(kR)\ A(k) - TlW(kR) A P (k) . (2.25) 
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The HI fluctuation is a sum of two parts, one which is correlated with the dark mat- 
ter distribution and an uncorrelated Poisson fluctuation Ap. The latter arises from 
the discrete nature of the HII regions and has a power spectrum n^ n . Also, W(y) = 
(3/y 3 )[siny — ycosy] is the spherical top hat window function arising from the Fourier 
transform of the spherical bubbles. This gives 

P A2 (k) = [1 _ Kb c W(kR)] 2 P{k) + i nW ( kR ^ 2 (2.26) 

and 

PA m (k) = [1 - Ub c W{kR)\ P{k) (2.27) 

which we use in equation (12. lip to calculate C;(Az/). Alternately, we have the HI power 
spectrum (Bharadwaj and Ali, 2005) 

P m (k) = 4i ([1 - Kb c W(kR) + f?P{k) + [W(fcjR)]2 | (2.28) 

[ ^HII J 

which we can use in equation (12.211) to calculate Ci(Ai>) in the flat-sky approximation. 

This model has three independent parameters, namely the average neutral fraction 
xm, the comoving radius of the ionized bubbles R and the bias of the bubble centers with 
respect to the dark matter b c . Our analysis assumes non-overlapping spheres and hence 
it is valid only when a small fraction of the HI is ionized and the bias is not very large. 
As a consequence we restrict these parameters to the range xhi > 0.5 and b c < 1.5. We 
note that in the early stages of re ionization (ie . Xhtt <S 1) equation (12.281) matches the 



HI power spectrum calculated by 



Wang fc Hul (120051 ) . though their method of arriving 



at the final result is somewhat different and is quite a bit more involved. 

Figure 12.11 shows the behaviour of Pni(k) for the two different models considered 
here. The cosmological parameters used throughout this paper are those determined as 
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Figure 2.1: The shows the HI power spectrum Pm(k) for the different models of the HI 
distribution considered here. Other than the one labeled DM, the curves are 
all for the PR model with different values of R and b c as indicated in the 
figure. Here the /i dependence has been incorporated by using the average 
value (/z 2 ) = 1/3. 



the best-fit values by W MAP 3-year data r elease, i.e., Q m = 0.23, VL^h 2 = 0.022 



0.96, h = 0.74, <t 8 = 0.76 (ISpereel et al. 



20061 ). Further, without any loss of generality, we 



have restricted our analysis to a single redshift z = 10 which corresponds to a frequency 
v = 129 MHz, and have assumed xm = 0.6 (implying 1Z — 2/3) which is consistent with 
currently favoured reionization models. 

The reason for choosing this particular redshift is that the effects of patchy reioniza- 
tion are expected to be most prominent around z ~ 10 in currently favoured reionization 
models. At higher redshifts, the reionization is in its preliminary stages (xhii 1) and 
the characteristic bubble size R is quite small. This implies that the effects of patchy 
reionization are not substantial and hence the HI distribution essentially traces the dark 
matter. Thus our results for the DM model are representative of what is expected at 
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higher redshifts, just that the overall normalization of Pni(k) would possibly be z depen- 
dent through the values of the growing mode of density perturbations and the neutral 
fraction. At lower redshifts the HI signal is expected to be drastically diminished be- 
cause most of the hydrogen would be ionized. Given this, it is optimum to study the 
HI signal properties at some intermediate redshift where xhii ~ 0.5 and P is reasonably 
large. For the currentl y favoured reionization scena rios, it seems that these properties 



2006al ). which we shall be studying in the 



are satisfied at z ~ 10 (jChoudhury fc Ferrara 
rest of this paper. 

The curve labeled DM in Figure 12.11 shows Pm{k) when the HI traces the dark 
matter. The characteristic scale in this power spectrum is set by the Fourier mode 
entering the horizon at the epoch of matter radiation equality. The imprint of the 
acoustic oscillations in the dark matter power spectrum holds interesting possibilities for 
determining cosmological parameters using high-z HI observations, we do not consider 
this here. The other curves in Figure 12.11 all show the PR model for different values of 
P and b c . The point to note is that for large values of the bubble size (P > 8 Mpc) 
the power spectrum is essentially determined by the Poisson fluctuation term Pni(k) ~ 
Xft{R 2 W 2 (kR)/n m i = xmiW 2 (kR)(47rR 3 /3), which scales as R 3 and is independent of 
the bias parameter b c . For length-scales larger than the bubble size (k < tt/R) we 
have W 2 (kR) ~ 1, and hence the power spectrum is practically constant Pni(k) ~ 
^hii(47T.R 3 /3). Around scales corresponding to the characteristic bubble size k « ir/R, 
the window function W(kR) starts decreasing which introduces a prominent drop in 
Pni(k). For smaller length-scales (k > tt/R), the power spectrum shows oscillations 
arising from the nature of the window function W(kR). At these scales, the amplitude 
of W 2 (kR) decreases as (/cP) -4 which is more rapid than P(k). Hence, at sufficiently 
large k the power spectrum Pni(k) is dominated by the dark matter fluctuations and 
it approaches the DM model, with the approach being faster for large P. We should 
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mention here that the oscillations in Pni(k) are a consequence of the fact that we have 
chosen the ionized bubbles to be spheres, all of the same size. In reality the ionized 
regions will have a spread in t he bubble shapes a nd sizes, and it is quite likely that such 



oscillation will be washed out (IWang fc Hu 



20051 ) but we expect the other features of the 



PR model discussed above to hold if the characteristic bubble size is large (R > 8 Mpc). 

For smaller values of R, the power spectrum Pni(k) could be dominated either by 
the term containing the dark matter power spectrum P(k) or by the Poisson fluctuation 
term, depending on the value of lZb c = 2b c /3. For values of b c « 3/2, the coefficient of 
P{k) in equation (I2.28f) tends to vanish, and hence the Poisson fluctuations dominate. 
This is obvious from the curve with parameters {R, b c } = {2 Mpc, 1.5} in Figure 12.11 
On the other hand, when the bias parameter is small, say b c < 1, the dark matter term 
dominates over the Poisson fluctuation term at large length-scales (k < 0.25 Mpc -1 ), 
which can be seen from the curve with {R,b c } = {2 Mpc, 1}. There is some difference 
between this curve and the DM model because of the TZb c W(kR) factor in equation 
(I2.28p . In fact, the the PR and DM models exactly coincide at large scales for {R, b c } = 
{2 Mpc, 0} which we have not shown separately in Figure 12.11 At small length-scales 
(k > 0.25 Mpc -1 ) the amplitude of the dark matter power spectrum becomes less than 
that of the Poisson fluctuation term, and hence Pni(k) is independent of b c (one can see 
that the curve having b c = 1 overlaps with the one having higher bias factor b c = 1.5). 
As mentioned earlier, at small length-scales (kR ^> 1) we expect the W 2 (kR) to decay 
rapidly as (&P)~ 4 , and as a consequence Pni(k) will basically trace the dark matter. 
It may be noted that for R = 2Mpc we do not notice this behaviour all the way till 
k = 10 Mpc -1 which is shown in the Figure 12.11 

In addition to the effects considered above, the random motions within clusters could 
significantly modify the signal by elongating the HI clustering pattern along the line of 
sight [the Finger of God (FoG) effect]. We have incorporated this effect by multiplying 
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Figure 2.2: The angular power spectrum of HI brightness temperature fluctuations for 
different models of the HI distribution assuming b c = 1 and a p = 0. 



the power spectrum P^Ck ) with an extra Lorentzian term (1 + k?,ap/a H ) (ISheth 



1996 



Ballinger et al 



19961 ) where o> is the one dimensional pair- velocity dispersion in 



relative galaxy velocities. 



2.3 Results 



We first consider the angular power spectrum Ci(Au) at Au = for which the results are 
shown in Figure [2~2l As discussed earlier, Ci(0) is essentially the 2D power spectrum of 
HI fluctuations evaluated at the 2D Fourier mode l/r v ~ I x 10 _4 Mpc _1 . The results for 
the DM model serve as the fiducial case against which we compare different possibilities 
for patchy reionization. 

For large bubble size {R > 8 Mpc) the HI signal is dominated by Poisson fluctuations 
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and it is well described through 



R 



+ l)Q oc v^iTT - I (2.29) 

r v 



on scales larger than the bubble. At these angular scales the HI signal is substantially 
enhanced compared to the DM model. For smaller bubble size, the large angle signal 
is sensitive to the bias b c . The signal is very similar to the DM model for b c = and 
it is suppressed for higher bias. In all cases (large or small bubble size), the signal is 
Poisson fluctuation dominated on scales comparable to the bubble size and it peaks at 
/ « irr u /R, with no dependence on b c . The HI signal traces the dark matter on scales 
which are much smaller than the bubble size. 

We next consider the behavior of Ki(Au), the frequency decorrelation function shown 
in Figure 12.31 For the DM model (upper left panel) where the HI fluctuations trace the 
dark matter we find that the frequency difference Au over which the HI signal remains 
correlated reduces monotonically with increasing I. For example, while for I = 100 
Ki(Au) falls to ~ 0.5 at Au ~ 500 KHz, it occurs much faster (Au ~ 10 KHz) for 
I = 10 5 . Beyond the first zero crossing Ki(Au) becomes negative (anti-correlation) and 
exhibits a few highly damped oscillations very close to zero. These oscillations arise 
from the cos term in equation ( 12. 21ft . The change in the behavior of ki(Au) for the DM 



Wang fc Hu 



(120051 ) 



model arising from the FoG effect is also shown in the same panel, 
have proposed that o v is expected to have a value ~ 30km/s at z ~ 8; in view of this, 
we show results for o p = 20and40km/s. We find that there is a discernible change at 
I > 10 4 , and the FoG effect causes the signal to remain correlated for a larger value of 
Au. For ij p = 20km/s, the change is at most 15% for I = 10 4 and around 100% at 
I = 10 5 . Though we have not shown it explicitly, we expect similar changes due to FoG 
effect in the PR model also. 
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Ay MHz Ay MHz 

Figure 2.3: The frequency decorrelation function ki(Au) [defined in equation (12.131) ] at 
I = 10 2 , 10 3 , 10 4 , 10 5 . Results are shown for the DM model and the PR model 
with b c = 1 and the R values shown in the figure. For the DM model, we 
show results incorporating the FoG effect using o> = 20and40km/s. For 
each / value K;(Az/) decreases faster for op = and slowest for ap = 40 km/s. 
There is a significant change due to the FoG effect only at / > 10 4 . 
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The patchy reionization model shows distinct departures from the DM model in the 
behavior of Ki(Av). This reflects the imprint of the bubble size and the bias on the 
Av dependence. For R = 2Mpc and b c = 1 (upper right panel) the large / (/ > 1000, 
comparable to bubble size) behavior is dominated by the Poisson fluctuation of the 
individual bubbles which makes ki(Av) quite distinct from the DM model. Notice that 
for I = 10 3 , Ki(Av) falls faster than the DM model whereas for I = 10 4 it falls slower than 
the DM model causing the I = 10 3 and 10 4 curves to nearly overlap. The oscillations 
seen in C\ as a function of I in Figure [2T21 are also seen in the Av dependence of K\{Av) 
at large I (10 5 ). The behavior at I = 10 2 is a combination of the dark matter and the 
ionized bubbles, and is sensitive to b c . For b c — 1, the initial decrease in Ki(Av) is 
much steeper than the DM model with a sudden break after which the curve flattens. 
Figure [231 shows the b c dependence for / = 10 2 and 10 3 . The bias dependence is weak for 
I = 10 3 where the Poisson fluctuations begin to dominate. For I = 10 2 , changing b c has a 
significant affect only near the break in Ki(Av) leaving much of the curve unaffected. For 
a smaller bubble size we expect a behavior similar to R = 2Mpc, with the b c dependence 
being somewhat more pronounced and the Poisson dominated regime starting from a 
larger value of I. 

For large bubble size (R > 8Mpc) the large angle HI signal (/ < nr^/R) is entirely 
determined by the Poisson fluctuations where the signal is independent of I. This is 
most clearly seen for R = 8Mpc where the Ki(Av) curves for I = 10 2 and I = 10 3 are 
identical. For both R = 8Mpc and 32Mpc the large / behavior of k,i(Au) approaches 
that of the DM model. 

In the final part we quantify the frequency difference Av across which the HI signal 
at two different frequencies remain correlated. To be more precise, we study the behavior 
of Av\ji which is defined such that Ki(Ai/u 2 ) = 1/2 ie. the correlation falls to 50% of 
its peak value at Av = 0. We study this for different angular scales (different I) for the 
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Figure 2.4: This shows the b c dependence of the frequency decorrelation function Ki(Au) 
for the PR model with R = 2 Mpc. 
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Figure 2.5: This shows Av\m vs. I for the DM model and the PR model with b c — 1 for 
the R values shown in the figure. 

various models of HI distribution considered here. The main aim of this exercise is to 
determine the frequency resolution that would be required to study the HI fluctuations 
on a given angular scale I. Optimally one would like to use a frequency resolution 
smaller than Az/i/2- A wider frequency channel would combine different uncorrelated 
signals whereby the signal would cancel out. Further, combining such signals would not 
lead to an improvement in the signal to noise ratio. Thus it would be fruitful to combine 
the signal at two different frequencies only as long as they are correlated and not beyond, 
and we use Avi/ 2 to estimate this. The plot of Aui/ 2 vs I for the different HI models is 
shown in Figure 12.51 

We find that for the DM model Avi/ 2 falls monotonically with I and the relation is 
well approximated by a power law 



Az/ 1/2 = 0.2 MHz x ( 



-0.7 



(2.30) 



34 



2.4 Implications tor Separating bignal from Foregrounds 

which essentially says that Avi/ 2 ~ 0.66 MHz on 1° angular scales, Av\j 2 ~ 0.04 MHz 
on 1' angular scales and Av\j 2 ~ 2 KHz on 1" angular scales. 

For the PR model, as discussed earlier, Ki(Av) is / independent on angular scales 
larger than the bubble size I < nr u /R. As a consequence Av\j 2 also is independent of I 
and it depends only on the bubble size R. This can be well approximated by 

Aui/2 ~ 0.04 MHz (2.31) 

which given a large value 1.3 MHz for R = 32Mpc while it falls below the DM model 
for 0.08 MHz for R = 2Mpc. The large I behavior of Av\j 2 approaches the DM model 
though there are oscillations which persist even at large I. 

We note that our findings are consistent with the earlier find ings of 



Bharad wai fc Ali 



(2005) whereas they significantly different from the results of ISantos. Cooray &: Knox 



(120051 ) who assume frequency channels of 1 MHz which is too large. 



2.4 Implications for Separating Signal from Foregrounds 



Astrophysical foregrounds are expected to be several order of magnitude stronger than 
the 21 cm signal. The MAPS foreg round contribution at a frequency v can be parametrized 



as (pantos. Cooray fc Knox 



2005|) 



{.Ay) 



(2.32) 



where Vf = 130 MHz and a is the mean spectral index. The actual spectral index varies 
with line of sight across the sky and this causes the foreground contribution to decorrelate 
with increasing frequency separation Au which is quantified through the foreground fre- 
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Figure 2.6: Angular power spectrum C/(0) at v = 129 MHz for the two most dominant 
foreground components, the diffuse galactic synchrotron radiation (GS) and 
the extragalactic point sources (PS) assuming S cut = 0.1 mJy. The expected 



signal is also shown for th e DM model and the PR. model with R 



quency decorrelation function Ii(Av) (jZaldarriaga. Furlanetto &: Hernquist 
has been modeled as 



6Mpc. 



2004 ) which 



Ii(Av) = exp 



- log? 1 



Au 



/2? 



(2.33) 



We consider the two most dominant foreground components namely extragalactic 
point sources and the diffuse synchrotron radiation from our own galaxy. Point sources 
above a flux level S cut can be identified in high-resolution images and removed. We as- 
sume S cu t = O.lmJy and adopt the parameter values from Table 1 of ' 



Santos. Cooray fc Knox 



(120051 ) for A, a, (3 and £. Figure |2~o1 shows the expected C;(0) for the signal and fore- 



grounds. The galactic synchrotron radiation dominates at large angular scales / < 10, 000 
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while the extragalactic point sources dominate at small angular scales. For all values of 
I, the foregrounds are at least two orders of magnitude larger than the signal. 

The foregrounds have a continuum spectra, and the contributions at a frequency 
separation Av are expected to be highly correlated. For Av = 1 MHz, the foreground 
decorrelation function Ii(Av) falls by only 2 x 10~ 6 for the galactic synchrotron radia- 
tion and by 3 x 10~ 5 for the point sources. In contrast, the HI decorrelation function 
ki(Av) is nearly constant at very small Av and then has a sharp drop well within 
1 MHz, and is largely uncorrelated beyond. This holds the promise of allowing the 
signal to be separated from the foregrounds. A possible strategy is to cross-correlate dif- 
ferent frequency channels of the full data which has both signal and foregrounds, and to 
use the distinctly different Av dependence to separate the signal from the foregrounds 



(jZaldarriaga. Furlanetto &: Hernquist 



20041 ) . An alternate approach is to subtract a 



20061 ) and then deter- 



best fit continuum spectra along each line of sight (iWang et all 
mine the power spectrum. This is expected to be an effective foreground subtraction 
method in data with very low noise levels. We consider the former approach here, and 
discuss the implications of our results. 

MAPS characterizes the joint I and Av dependence which is expected to be different 
for the signal and the foregrounds. For a fixed I, it will be possible to separate the 
two with relative ease at a frequency separation Av if the decrement in the signal 
Ci(0)[l — Ki(Av)] is more than that of the foregrounds Cj(0)[l — Ii(Av)]. Note that 
because the foregrounds are much stronger than the HI signal, a very small decorrelation 
of the foreground contribution may cause a decrement in Ci(Av) which is larger than 
that due to the signal. We use Fi(Av) defined as the ratio of the two decrements 

F (Ay) = {a(0)[l-^(A^)]}si g nal , , 

grounds 
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Figure 2.7: This shows F\(Av) (defined in equation 12.341) for the Av values shown in the 
figure. We consider both the DM model(left panel) and the PR model(right 
panel). 

to asses the feasibility of separating the HI signal from the foregrounds. This gives an 
estimate of the accuracy at which the Av dependence of the foreground Ci(Av) has to be 
characterized for the signal to be detected. Note we assume that the ( u + f Al/ ) a term in eq. 
( I2.32p can be factored out before considering the decrement in the foreground. Figure 
12.71 shows the results for the DM model and the PR model with R — 16 Mpc. First 
we note that Fi(Av) peaks at the angular scales corresponding to I ~ 10,000 (ie. 2') 
and the prospects of separating the signal from the foregrounds are most favorable 
at these scales. A detection will be possible in the range I > 1000, Av < lOKHz 
and / > 400, Av < lOOKHz for the DM and PR models respectively provided the Av 
dependence of the foregrounds Ci(Av) can be characterized with an uncertainty less than 
order unity. The I and Av range would increase if the Av dependence of the foreground 
C[(Av) were characterized to 10% accuracy. The largest angular scales ( / < 100 ) would 
require an accuracy better than 1% which would possibly set the limit for forthcoming 
observations. 

The angular modes / = 1,000 and / = 10,000 correspond to baselines with antenna 
separations of ~ 300 m and ~ 3 km respectively. This baseline range is quite well 
covered by the GMRT, and also the forthcoming interferometric arrays. This is possibly 
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the optimal range for a detection. A possible detection strategy would be to use the 
Av behavior of Ci(Av) in the range where Fi(Av) <C 1 to characterize the foreground 
contribution. This can be extrapolated to predict the foreground contribution at small 
Av and any excess relative to this prediction can be interpreted as the HI signal. A very 
precise determination of the Av dependence of the foreground contribution would require 
a very large Av range in the region where Fi(Av) <C 1, and a bandwidth of ~ 10MHz 
would be appropriate. On the other hand, at I ~ 10,000 the HI Ci(Av) decorrelates 
within ~ 50KHz [or equivalently Fi(Av) shows a considerable drop between 10 KHz and 
100 KHz, see Figure E7j, and it would be desirable to have a frequency resolution better 
than ~ lOKHz to optimally differentiate between the signal and the foregrounds. A 
lower resolution of ~ 20 KHz would possibly suffice at I ~ 1, 000, particularly if the PR 
model holds. 
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3 Detecting Ionized Bubbles in 



Redshifted 21 cm Maps 



3.1 Introduction 

In this Chapter we consider the possibility of detecting ionized bubbles in redshifted 
21 cm HI maps. An ionized bubble embedded in HI will appear as a decrement in the 
background redshifted 21 cm radiation. This decrement will typically span across several 
pixels and frequency channels in redshifted 21 cm maps. Detecting this is a big chal- 
lenge because the HI signal (~ 1 mJy or lower ) will be buried in foregrounds which are 
expected to be at least 2 — 3 orders of magnitude larger. An objective detection criteria 
which optimally combines the entire signal in the bubble while minimizing contributions 
from foregrounds, system noise and other such sources is needed to search for ion- 
ized bubbles. The noise in diffe rent pixels of maps obtained fr o m rad io-interferometric 



observations is correlated (eg. 



Thompson. Moran Swensonl (119861 )). and it is most 



convenient to deal with visibilities instead. These are the primary quantities that are 
measured in radio-interferometry. 



This chap ter is adapted from the paper "D e tectin g ionized bubbles in redshifted 21 cm 
maps" bv \Datta. Bharadwai & Choudhuru (200'a) . 
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In this Chapter we develop a visibility based formalism to detect an ionized bubble 
or conclusively rule it out in radio-interferometric observations of HI at high redshifts. 
We apply our formalism for detecting ionized bubbles to make predictions for the GMRT 
and for one of the forthcoming instruments, namely the MWA. For both telescopes we 
investigate the feasibility of detecting the bubbles, and in situations where a detection 
is feasible we predict the required observation time. For both telescopes we make pre- 
dictions for observations only at a single frequency (150MHz), the aim here being to 
demonstrate the utility of our formalism and not present an exhaustive analysis of the 
feasibility of detecting ionized bubbles in different scenarios and circumstances. For the 
GMRT we have used the telescope p arameters from t heir w ebsite, while for the MWA 
we use the telescope parameters from 



Bowman et al. 



(120061 ) 



The outline of the Chapter is as follows: In Section 3.2 we discuss various sources 
which are expected to contribute in low frequency radio-interferometric observation, 
this includes the signal expected from an ionized bubble. In Section 3.3 we present the 
formalism for detecting an ionized bubble, and in Section 3.4 we present the results and 
discuss its implications. The cosmological parameters used throughout this Chapter are 
those determined as the best-fit mean values for a flat ACDM model by WMAP 3-year 



data release, i.e., Q m = 0.23, Vl b h 2 = 0.022, n s = 0.96, h = 0.74, a 8 = 0.76 flSpergel et al 
20061 ) . 



3.2 Different Sources that Contribute to Low Frequency 
Radio Observations 

The quantity measured in radio-interferometric observations is the visibility V(U, v) 
which is measured in a number of frequency channels v across a frequency bandwidth 
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B for every pair of antennas in the array. For an antenna pair, it is convenient to use 
U = d/X to quantify the antenna separation d projected in the plane perpendicular to 
the line of sight in units of the observing wavelength A. We refer to U as a baseline. 
The visibility is related to the specific intensity pattern on the sky I v (9) as 



V(U, v) = J d 2 BA{B)I v 0)e 2M (3.1) 



where 9 is a two dimensional vector in the plane of the sky with origin at the center of 
the field of view, and A(8) is the beam pattern of the individual antenna. For the GMRT 
this can be well approximated by Gaussian A{6) = e~ d2 ^ e ° 2 where 6$ ~ 0.6 #fwhm and 
we use the values 2.28° for O at 150 MHz for the GMRT. Each MWA antenna element 
consists of 16 crossed dipoles distributed uniformly in a square shaped tile, and this is 
stationary with respect to the earth. The MWA beam p attern is quite compli cated, 



2003). Our 



and it depends on the pointing angle relative to the zenith (IBowman et al. 
analysis largely deals with the beam pattern within 1° of the pointing angle where it is 
reasonable to approx imate the beam as being circularly symmetric (Figures 3 and 5 of 



Bowman et al. 



20071 ). We approximate the MWA antenna beam pattern as a Gaussian 
with 9q = 18° at 153 MHz. Note that the MWA primary beam pattern is better modeled 
as A{9) oc cos 2 (K6), but a Gaussian gives a reasonable approximation in the center of 
the beam which is the region of interest here. Equation (13.1 ft is valid only under the 
assumption that the field of view is small so that it can be well approximated by a plane, 
or under the unlikely circumstances that all the antennas are coplanar. 

The visibility recorded in 150 MHz radio-interferometric observations is a combina- 
tion of three separate contributions 

V(U, v) = S(U, v) + N(U, v) + F{U, v) (3.2) 
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— * — * 

where S(U, v) is the HI signal that we are interested in, N(U, v) is the system noise which 
is inherent to the measurement and F(U, v) is the contribution from other astrophysical 
sources referred to as the foregrounds. Man-made radio frequency interference (RFI) 
from cell phones and other communication devices are also expected to contribute to 
the measured visibilities. Given the lack of a detailed model for the RFI contribution, 
and anticipating that it may be possible to remove it before the analysis, we do not take 
it into account here. 

3.2.1 The HI signal from ionized bubbles 

According to models of reionization by UV sources, the early stages of reionization are 
characterized by ionized HII regions around individual source (QSOs or galaxies). As a 
first approximation, we consider these regions as ionized spherical bubbles characterized 
by three parameters, namely, its comoving radius Rb, the redshift of its center z c and the 

— * 

position of the center determined by the two-dimensional vector in the sky-plane 9 C . The 
bubble is assumed to be embedded in an uniform intergalactic medium (IGM) with a 
neutral hydrogen fraction xm- We use r v to denote the comoving distance to the redshift 
where the HI emission, received at a frequency v = 1420 MHz/ (1 + z), originated, and 
define r' u = dr v jdv. The planar section through the bubble at a comoving distance 
r v is a disk of comoving radius R u = Rb^l — (Au/ Aub) 2 where Av = v c — v is the 
distance from the the bubble center v c in frequency space with v c = 1420MHz/(l + z c ) 
and Aub = Rb/r' Uc is the bubble size in the frequency space. The bubble, obviously, 
extends from v c — Avb to v c + Avb in frequency and in each frequency channel within 
this frequency range the image of the ionized bubble is a circular disk of angular radius 
6 V = R u /r u ; the bubble is not seen in HI beyond this frequency range. Under such 
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assumptions, the specific intensity of the redshifted HI emission is 



I v {9) = I v x m 



1 



e-e c 



ei 



V — V r 



(3.3) 



where l u = 2.5 x 10 2 ^ {jq^j (tt) \W(z)) * s ^ ne ra diation background from the uniform 
HI distribution and 0(x) is the Heaviside step function. 

The soft X-ray emission from the quasar responsible for the ionized region is expected 



to heat the neutral IGM in a shell around the ionized bubble. The H 



this shell is expected to be somewhat higher than l v ( iWyithe &: Loeb . 



emission from 



200481 ). We do 



not expect this to make a very big contribution, and we do not consider this here. 

If we assume that the angular extent of the ionized bubble is small compared to the 
angular scale of primary beam ie. 9 V <C 9q, we can take A(9) outside the integral in eq. 
(13. ip and write the signal as A(6 C ) J d 2 9I u (9)e 2me ' u , which essentially involves a Fourier 
transform of the circular aperture ^1— | 9 — 9 C \ r v jR^. For example, a bubble of 
radius as large as 40 Mpc at z = 8.5 would have an angular size of only Q v « 0.25° which 
satisfies the condition 6 U <C 9 . In a situation where the bubble is at the center of the 
field of view, the visibility is found to be 



S centeT (U,v) = -nl u x m 0l 



2J 1 {2nU6 v ) 
2nU6„ 



1 



V — V r 



(3.4) 



where J\{x) is the first order Bessel function. Note that Center (U, v) is real and it is the 
Fourier transform of a circular aperture. The uniform HI background also contributes 
7 l/ 7r^o6 _7r26 '° !72 to the visibility, but this has been dropped as it is quite insignificant at 
the baselines of interest. Note that the approximations used in eqs. ( 13.41) have been 
tested extensively by comparing the values with the numerical evaluation of the integral 
in eq. ( 13. If) . We find that the two match to a high level of accuracy for the situations of 
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Figure 3.1: Signal from a spherical ionized bubble of comoving radius lOMpc as a func- 
tion of baseline U for different frequency channels. 

interest here. In the general situation where the bubble is shifted by 6 C from the center 
of the field of view, the visibility is given by 



S(U, v) = e- e2 ^oe 2 ^S centei (U, u) 



(3.5) 



i.e., there is a phase shift of e 27nC/ ' 0c and a e 6c ^ e ° drop in the overall amplitude. 

Figures 13.11 and 13.21 show the U and Av dependence of the visibility signal from an 
ionized bubble with Rb = lOMpc located at the center of the field of view at v c = 
153 MHz (z c = 8.3), assuming xhi = 1. The signal extends over Au = ±Az/& where 
Ai>b = 0.56 MHz. The extent in frequency Av\, = Rb/r' Uc scales oc Rb when the bubble 
size is varied. The Bessel function J\{x) has the first zero crossing at x — 3.83. As a 



result, the signal S(U,u) extends to U = 0.61r y [_R b A/l — (Av/Avb) 2 ] 1 where it has 
the first zero crossing, and Uq scales with the bubble size as Uq oc 1/Rb- The peak value 
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Figure 3.2: Signal from a spherical ionized bubble of comoving radius 10 Mpc as a func- 
tion of Av = v — v c for different baselines. 
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of the signal is 5(0, v) = -Kl v (Rb/r u ) 2 yJ\ — (Av/Ais b ) 2 and scales as 5(0, v) oc R\ if the 
bubble size is varied. We see that the peak value of the signal is 5(0, is c ) = 70 /xJy for 
bubble size Rb = lOMpc and would increase to 1.75 mJy if Rf, = 50Mpc. Detecting 
these ionized bubbles will be a big challenge because the signal is buried in noise and 
foregrounds which are both considerably larger in amplitude. Whether we are able to 
detect the ionized bubbles or not depends critically on our ability to construct optimal 
filters which discriminate the signal from other contributions. 



3.2.2 HI fluctuations 

In the previous sub-section, we assumed the ionized bubble to be embedded in a per- 
fectly uniform IGM. In reality, however, there would be fluctuations in the HI distri- 
bution in the IGM which, in turn, would contribute to the visibilities. This contri- 
bution to the HI signal can be treated as a random variable S(U, v) with zero mean 
(S(U, u)) = 0, whose statistical properties are characterized by the two-visibility corre- 



lation (S(Ui,vi)S(U2,U2)). This is r elated to Pm(k) t 



radiation efficiency in redshift space (IBharadwaj &: Ali 



r e pow er spectrum of the 21 cm 



20041 ) through 



2r 2 

V 



POO 

: / dfc|| PHi(k) cos^nr^Az/) 
Jo 



(3.6) 



where 5jj jj 2 is the Kronecker delta it. different baselines are uncorrelated, To estimate 
the contribution from the HI fluctuations we make the simplifying assumption that the 
HI traces the dark matter, which gives Pm(k) = xfjj (1 + u 2 ) 2 P(k) where P(k) is the 
dark matter power spectrum and \i is the cosine of the angle between k and the line of 
sight. This assumption is reasonable because the scales of interest are much larger than 
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the Jeans length Aj ~ 10 — 100 kpc, and we expect the HI to cluster in the same way as 
the dark matter. 

In addition to the above, there could be other contributions to the HI signal too. 
For example, there would be several other ionized regions in the field of view other 
than the bubble under consideration. The Poisson noise from these ionized patches will 
increase the HI fluctuations and there will also be an overall drop in the contribution 
because of the reduced neutral fraction. These effects will depend on the reionization 
model, and the simple assumptions made would only provide a representative estimate 
of the actual contribution. Figure 13.31 shows the expected contribution from the HI 
fluctuations (HF) to the individual visibilities for GMRT and MWA. Note that while 
this can be considerably larger than the signal that we are trying to detect (particularly 
when the bubble size is small), there is a big difference between the two. The signal 
from the bubble is correlated across different baselines and frequency channels whereas 
the contribution from random HI fluctuations is uncorrelated at different baselines and 



it become uncorrelated beyond a certain fre quency separation Az/ (IBharadwaj fc Ali 



2005 : 



Datta. Choudhury &: Bharadwaj 



20071 ) 



3.2.3 Noise and foregrounds 



The system noise contribution N(U, v) in each baseline and frequency channel is ex- 
pected to be an independent Gaussian random variable with zero mean ((N) = 0) and 
whose variance is independent of U a nd u c . The predicted rms. noise contribution is 



([Thompson. Moran fc Swensonl (119861 1) 



(N 2 ) 



V2k B T, 



sys 



(3.7) 
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where T sys is the total system temperature, k B is the Boltzmann constant, A e ff is the 
effective collecting area of each antenna, Av c is the channel width and At is the correlator 
integration time. Equation ( 13. 7ft can be rewritten as 

-1/2 / A , \ -1/2 



where C x varies for different interferometric arrays. Using the GMRT parameters T sys 
482K and A eff /2k B = 0.33 K/Jy at 153MHz gives C x = 1 .03Jy for the GMRT where 



200J) gives 



as for MWA T sys = 470K and A eff /2k B = 5 x 10~ 3 K/Jy (iBowman et al. 
C x = 65.52Jy. The rms noise is reduced by a factor a/ At/t b s if we average over t b s / At 
independent observations where t f, s is the total observation time. Figure 13.31 shows the 
expected noise for a single baseline at 153MHz for Av c = 50 KHz and an observation 
time of 100 hrs for both the GMRT and MWA. Though T sys is nearly equal for the GMRT 
and the MWA, the noise in a single baseline is expected to be 60 times larger for MWA 
than that for the GMRT. This is a because the individual antennas have a much larger 
collecting area at the GMRT as compared to the MWA. The fact that the MWA has 
many more antennas (N = 500) as compared to the GMRT (N = 30) compensates for 
this. Note that nearly half (16) of the GMRT antennas are at very large baselines which 
are not particularly sensitive to the signal on the angular scales the ionized bubble, and 
only the other 14 antennas in the 1 km x 1 km central square will contribute towards 
detecting the signal. For both the GMRT and the MWA, T sys is dominated by the sky 
contribution T s k y with the major contribution coming from our Galaxy. We expect T sys 
to vary depending on whether the source is in the Galactic plane or away from it. The 
value which we have used is typical for directions off the Galactic plane. Further, the 
noise contribution will also be baseline dependent which is not included in our analysis. 
Contributions from astrophysical foregrounds are expected to be several order of 
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magnitude stronger than the HI signal. Extragalactic point sources and synchrotron 
radiation from our Galaxy are predicted to be the most dominant foreground compo- 
nents. Assuming that the foregrounds are randomly distributed, with possible cluster- 
ing, we have (F(U,v)) = for all the baselines other than the one at zero spacing 
(U = 0), which is not considered in this work. The statistical properties are charac- 
terized by the two-visibility correlation (F(Ux, vi)F(U2, Vz))- We express this (details 
in Appendix E2) in terms of the multi-frequency angular power spectrum (hereafter 
MA PS) C][v-\ v?) of the brightness temperature fluctuations at the frequen cies v\ and v 2 



as (jSantos. Cooray &: Knox 



2005 



Datta. Choudhury &: Bharadwaj 



20071 ) 



m 
n + 

( 9 4) ( 9 4) <W**)- (3.9) 



\&r)„ \df) 



where (dB/dT) u = 2kBV 2 /c 2 is the conversion factor to specific intensity, and we have 
assumed that the primary beam pattern A{6) = e~ e2 ^ 9 ° is frequency dependent through 
8o oc v~ x and use Q\ and 82 to denote the value of #0 at v\ and v 2 respectively. Note that 
the foreground contribution to different baselines are expected to be uncorrelated. 
For each component of the foreground the MAPS is modeled as 

Cl(viV2 ) = A {^f (^-J (3.10) 

where Vf = 130 MHz, and for each foreground component A, (3 and a are the ampli- 
tude, the power law index of the angular power spectrum and the mean spectral index 
respectively. The actual spectral index varies with line of sight across the sky and this 
causes the foreground contribution to decorrelate with increasing frequency separation 
IS.V = \vi — which is quantified through the foreground frequency decorrelation func- 
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tion z/ 2 ) (jZaldarriaga. Furlanetto & Hernquist 



20041 ) which has been modeled as 



h(y\ V2) = exp 



-log? (-)/2£ 2 



(3.11) 



We consider the two most dominant foreground components namely extragalactic point 
sources and the diffuse synchrotron radiation from our own galaxy. Point sources above 
a flux level S cu t can be identified in high-resolution continuum images and removed. We 
note that absence of large baselines at the MWA restricts the angular resolution, but 
it may be possible to use the large frequency bandwidth 32 MHz to identify continuum 
point sources in the frequency domain. S cut depends on a the rms. noise in the im- 
age. We use S cu t = 5cr where a is the rms noise in the image given by (assuming 2 
polarizations) 



a 



B \ 



-1/2 



^/2N~ b ViMHzy V lsec , 
where iV b = iV(iV — l)/2 is the number of independent baselines, N is the number of 
antennas in the array, B is the total frequency bandwidth and t b s the total observation 
time. For t D b s = 100 hrs and B = 6 MHz we have S cu t = 0.1m Jy for the GMRT and using 
B = 32 MHz it gives S cut = 0.2 mJy for the MWA. The value of S cut will be smaller for 
longer observations, but reducing S cu t any further does not make any difference to our 
results so we hold S cu t fixed at these values for the rest of our analysis. The confusion 
noise from the unresolved point sources is a combination of two parts, the Poisson 
contribution due to the discrete nature of these sources and the clustering contribution. 
The amplitude of these two contributions have different S cut dependence. The parameter 



tobs 



-1/2 



(3.12) 



values that we have used are listed in 
values from 



DiMatteo. et al. 



'able 13.11 We have adopted the parameter 



Santos. Coorav fc Knox! (120051 ) and incorporated the S cut dependence from 



200^ 
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Table 3.1: Parameters values used for characterizing different foreground contributions 



Foregrounds 


A(mK 2 ) 


a 


P 




Galactic synchrotron 


700 


2.80 


2.4 


4 


Point source 
(clustered part) 


61 (o£3y) 


0.5 


2.07 


1.1 


2 


Point source 
(Poisson part) 


°-16 (o.tly, 


1.25 


2.07 





1 



Figure 13.31 shows the expected foreground contributions for the GMRT and MWA. 
The galactic synchrotron radiation is the most dominant foreground component at large 
angular scales (U < 1000 for GMRT and U < 2000 for MWA), while the clustering of 
the unresolved extragalactic point sources dominates at small angular scales. For all 
values of U, the foregrounds are at least four orders of magnitude larger than the signal, 
and also considerably larger than the noise. 

The MWA has been designed with the detection of the statistical HI fluctuation signal 
in mind, and hence it is planned to have a very large field of view. The foreground 
contribution to a single baseline is expected to be 10 times stronger for the MWA 
than for the GMRT because of a larger field of view. As we shall show later, the 
increased foreground contribution is not a limitation for detecting HII bubbles. The 
foregrounds have a continuum spectra, and the contribution at two different frequencies 
at a separation Au are expected to be highly correlated. For Au = 1 MHz, the foreground 
decorrelation function Ii(Au) falls by only 2 x 10 -6 for the galactic synchrotron radiation 
and by 3 x 10 -5 for the point sources. In contrast, the signal from an ionized bubble peaks 
at a frequency corresponding to the bubble center and falls rapidly with Au (Figure |3~2"1) . 
This holds the promise of allowing the signal to be separated from the foregrounds. 
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Figure 3.3: The magnitude of the different contributions to the visibility V(U, v) at v — 
153 MHz as a function of U. The signal, foregrounds (FG), noise (NS) and 
HI fluctuations (HF) contributions are shown for the GMRT (left) and MWA 
(right). The expected signal is shown for bubbles with radius R = 10 Mpc 
and R = 50 Mpc. The noise is estimated for a single baseline assuming an 
observation time t b s = 100 hrs and channel width Au c = 50 KHz. 

3.3 Formalism for Detecting the Ionized Bubble 



We consider a radio- interferometric observation of duration t f, s , carried out over the 
frequency range v\ to v-i- The HI signal from an ionized bubble, if it is present in the 
data, will be buried in foregrounds and noise both of which are expected to be much 
larger. In this Section we present a filtering technique aimed at detecting the signal 
from an ionized bubble if it is present in our observations. To detect the signal from an 
ionized bubble of radius Rf, with center at redshift z c (or frequency v c ) and at an angle 
9 C from the center of the field of view, we introduce an estimator E[Rb, z c , 9 C ] defined as 



E 



Y,S}{U a ,v b )V{U a ,v b ) 

. ct,b 



I 



a.b 



(3.13) 



where Sf(U,i/) is a filter which has been constructed to detect the particular ionized 
bubble. Here U a and u b refer to the different baselines and frequency channels in our 
observations, and in eq. (I3.13P we are to sum over all independent data points (visibili- 
ties). Note that the estimator E and the filter Sf(U, v) both depend on [Rj,, z c , 9 C ], the 
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parameters of the bubble we wish to detect, but we do not show this explicitly. The 
values of these parameters will be clear from the context. 

We shall be working in the continuum limit where the two sums in eq. (13 . 1 3|) can be 
replaced by integrals and we have 

E= J d 2 U J dup N (U,u) S f *{U,v)V{U,v) (3.14) 

d 2 Udv Pn{U , v) is the fraction of data points ie. baselines and frequency channels in the 
interval d 2 U dv. Note that Pn{U, v) is usually frequency dependent, and it is normalized 
so that J d 2 U J dv p N {U ,u) = 1. We refer to Pn(U,i>) as the normalized baseline 
distribution function. 

We now calculate (E) the expectation value of the estimator. Here the angular 
brackets denote an average with respect different realizations of the HI fluctuations, 
noise and foregrounds, all of which have been assumed to be random variables with zero 
mean. This gives (V(U, is)) = S(U, is) and 

(E) = J d 2 U J dvp N (U,v) S f *(U,v)S{U,v) (3.15) 

We next calculate the variance of the estimator which is the sum of the contributions 
from the noise (NS), the foregrounds(FG) and the HI fluctuations (HF) 

<(A£) 2 > = ((E — (E)) 2 } 



(AE) } +((AE) 2 ) +((AEy 

I NS \ / FG \ / HF 

(3.16) 

To calculate the noise contribution we go back to eq. (13.131) and use the fact that the 
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noise in different baselines and frequency channels are uncorrected. We have 



((AE) 2 ) m = (N 2 ) 



I S f (U a} is b ) 



a, b 



(3.17) 



which in the continuum limit is 



((AE) 2 ) 



NS 



<iv 2 >/£i 



o,6 



x / d 2 U / dvp N (U,v) | 3 f (U,u) 



(3.18) 



The term 
then have 



(N 2 )/ z : 2 ab 1 is the same as ex, the rms. noise in the image (eq. I3.12p . We 
((AE) 2 ) NS = a 2 [ d 2 U [ dvp N {U,v) | S f {U,v) | 2 . (3.19) 



For the foreground contribution we have 



(AE) 2 ) 



FG 



d 2 Ux d 2 U 2 / di/i / dz/ 2 



xpjv(C/i, ^i)PTv(f/ 2 , i^)S>*(Z7 1} vx)S f *{U % , u 2 ) 
In the continuum limit we have (details given in Appendix [C]) 



z/ 2 



(3.20) 



(3.21) 
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which gives the variance of the foreground contribution to be 



xp N (U, u l )p N (U, u 2 )S f *(U, u 1 )S f (U, u 2 ) 



xC w (^,!/ 2 ) (3.22) 



We use eq. (13.221) to calculate ( (AE) 2 ) too, with the difference that we use the 

\ / HF 



power spectrum C 2 ^u (u, z/+Az/) for the HI fluctuation from 



Datta. Choudhury fc Bharadwaj 



(120071 ) instead of the foreground contribution. 



In an observation it will be possible to detect the presence of an ionized bubble having 



parameters [Rj,, z c , 8 C ] at, say 3-sigma confidence level, if (E) > 3y ((AE) 2 ). In such 
a situation, an observed value E can be interpreted as a detection with 99.7% (i.e., 



3-sigma) confidence if E Q > 3y ((AE) 2 ). The presence of the ionized bubble can be 



ruled out at the same level of confidence if (E) — E D > 3 y ((AE) 2 ). 



3.3.1 Baseline distribution 

In this subsection we discuss the normalized baseline distribution function pn(U,v) 
which has been introduced earlier. Figure 13.41 shows the baseline coverage for 14 hrs 
of observation towards a region at declination 5 = 45° with the GMRT at 153MHz. 
In this figure u and v refer to the Cartesian components of the baselines U. Note 
that the baseline distribution is not exactly circularly symmetric. This asymmetry 
depends on the source declination which would be different for every observation. We 
make the simplifying assumption that the baseline distribution is circularly symmetric 
whereby Pn(U, v) is a function of U. This considerably simplifies our analysis and gives 
reasonable estimates of what we would expect over a range of declinations. Figure 13.51 
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Figure 3.4: This shows the baseline coverage for 14hrs of GMRT 153 MHz observation 
at 45° declination. 
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Figure 3.5: This shows the normalized baseline distribution Pn{U, v) for the GMRT and 
the MWA at 153 MHz. The wiggly curve shows the actual values for the 
GMRT observation shown in Figure [33 and the smooth curve is the analytic 
fit. 
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shows pn{U,u) for the GMRT determined from the baseline coverage shown in Figure 
I3.4L We find that this is well described by the sum of a Gaussian and an expon ential 



2003 ) with 



distribution. The GMRT has a hybrid antenna distribution (jChengalur et al.l . 
14 antennas being randomly distributed in a central square approximately 1 km x 1 km 
and 16 antennas being distributed along a Y each of whose arms is 14 km long. The 
Gaussian gives a good fit at small baselines in the central square and the exponential 
fits the large baselines. Determining the best fit parameters using a least square gives 



Pn(U, v) 



1 

B 



A 



1 km^ 
9.70 x 10~ 3 exp 



0.21 exp 



U 2 X 2 



2a 2 



UX-b 
cT~ 



(3.23) 



where a = 0.382 km, b = 0.986 km, d = 3.07km and B is the frequency bandwidth which 
has a maximum value of 6 MHz. 
Following 



Bowman et al. 



(120061 ) we assume that the MWA antennas are distributed 
within a radius of 0.750 km with the density of antennas decreasing with radius r as 
Pant{ r ) oc r~ 2 and with a maximum density of one antenna per 18m 2 . The normalized 
baseline distribution is estimated in terms of p an t(r) and we have 



Pn{U,v) 



X 



x 



4.4 x 10 2 B VI km 
/ Pant(\r- A U Qd0 

J<h=Q 



d 2 rp ant (r) 



r=0 



(3.24) 



where the bandwidth B is 32 MHz, \r - A U\ = (r 2 + U 2 A 2 - 2r A Ucos<fi) 1/2 . Note that 
Pn(U, v) depends on the observed frequency. Figure 13751 shows the normalized baseline 
distribution function Pn{U, u) for both the GMRT and the MWA. We see that maximum 
baseline for the GMRT is U max ~ 10, 000 whereas U max ~ 750 for the MWA. However, 
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the smaller baselines will be sampled more densely in the MWA as compared to the 
GMRT. 



3.3.2 Filter 

It is a major challenge to detect the signal which is expected to be buried in noise and 
foregrounds both of which are much stronger (Figure I3~3"j) . It would be relatively simple 
to detect the signal in a situation where there is only noise and no foregrounds. The 
signal to noise ratio (SNR) is maximum if we use the signal that we wish to detect as 
the filter (ie. Sf(U, v) = S(U, v)) and the SNR has a value 



((E)) 



((A£) 2 ) ns 



a 



~l 0.5 



d 2 U / dvp N (U,v) | S(U,u) 



(3.25) 



The observing time necessary for a 3-a detection (i.e., SNR = 3) would be the least 
for this filter. Note that the factor —rrxmlv outside the signal (eq. 13.41) is almost 
constant along the line of sight of the ionized bubble. This factor does not affect the 
value of the quantity SNR= (( E J> =. One can drop this term from the filter Sf(U, v) 
without losing the effectiveness of the method. The difficulty with using this filter is 
that the foreground contribution to J ((AE) 2 ) is orders of magnitude more than ((E)). 
The foregrounds, unlike the HI signal, are all expected to have a smooth frequency 
dependence and one requires filters which incorporate this fact so as to reduce the 
foreground contribution. We consider two different filters which reduce the foreground 
contribution, but it occurs at the expense of reducing the SNR, and t Q b s would be more 
than that predicted by eq. (13.251) . 

The first filter (Filter I) subtracts out any frequency independent component from 
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the frequency range u c — B' /2 to u c + B /2 with B' < B ie. 



S f (U, u) 



A, 



S(U, u) 



9(1 - 2 | v - v c | /£') r c+B '/2 



B' 



S{U, v') dv' 



u c -B /2 



(3.26) 



This filter has the advantage that it does not require any prior knowledge about the 
foregrounds except that they have a continuous spectrum. It has the drawback that 
there will be contributions from the residual foregrounds as all the foregrounds are 
expected to have a power law spectral dependence and not a constant. A larger value of 
B' causes the SNR to increases, and in the limit B — » 00 the SNR approaches the value 
given in eq. (13.251) . Unfortunately the residues in the foregrounds also increase with B' . 
We use B' = 4A^ provided it is less than B, and B' = B otherwise. 

The frequency dependence of the total foreground contribution can be expanded in 
Taylor series. Retaining terms only up to the first order we have 



Ci{pi, v 2 ) = Ci(u c , v c ) [1 - (Az/i + Au 2 ) a eff /is c 



(3.27) 



where Au = v — v c and a e ff = , A , < " 1000 /!l. * is the effective spectral index, here i refers 
to the different foreground components. Note that a e // is / dependent. The second filter 
that we consider (Filter II) allows for a linear frequency dependence of the foregrounds 
and we have 



S f (U, v)=(l + a eff Av/v c ) 



S(U, v) 



9(1-2 I v-v c \ fB') 
B 



Vc+B /2 



v c -B' /2 



S(U, v) di/ 
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(3.28) 

Note that for both the filters we include an extra factor (A c /A) 2 . This is introduced with 
the purpose of canceling out the A 2 dependence of of the normalized baseline distribution 
function Pn(U, v) and this substantially reduces the foreground contribution. 



3.4 Results and Discussions 



We first consider the most optimistic situation where the bubble is at the center of 
the field of view and the filter center is exactly matched with the bubble center. The 
size distribution of HII regions are quite uncertain, and would depend on the reion- 
ization history and the distribution of ionizing sources. However, there are some in- 
dicatio ns in the literature on what co uld be the typical size of HII regions. For ex- 



ample, 



Wyithe. Loeb fc Barnes! (120051 ) deduce from proximity zone effects that Rb 



35Mpc at z ~ 6, which should be con sidered as a lower limit. On the other hand, 



Furlanetto. McQuinn fc Hernquist 



(120061 ) (Figure 1(a)) infer that the characteristic bub- 



ble size Rb > 10 Mpc at z — 8 if the ionized fraction Xi > .75 (Rh ~ 50Mpc if Xj y 0.9 



20071 ) 



Theoretical models which match a variety of observations (IChoudhury fc Ferrara 
imply that could be as high as 90 % at z ~ 8, which would mean bubble sizes of 
~ 40 — 50 Mpc. To allow for the large variety of possibilities, we have presented results 
for a wide range of Rb values from 2 Mpc to 50 Mpc. We restrict our analysis to a 
situation where the IGM outside the bubble is completely neutral (xhi = 1)- The sig- 
nal would fall proportional to xm if the IGM outside the bubble were partially ionized 



(xhi < 1). The expected signal (E) and 3-sigma fluctuation 3 x y ((AE) 2 ) from each of 
the different components discussed in Sections 2 and 3 as a function of bubble size Rb 
are shown in Figures l3~6l and I3TT1 Both the figures show exactly the same quantities, the 
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Figure 3.6: The signal quantified through the expectation value of the estimator (E) 
for Filter I. The other components (NS - Noise, FG - Foregrounds, HF - 
HI Fluctuations) are quantified through their contribution to the 3-sigma 



fluctuation 3 x y ((AE) 2 ). 
only difference being that they refer to Filter I and Filter II respectively. A detection is 



possible only in situations where (E) > 3 x y ((AE) 2 }, the rhs. now refers to the total 
contribution to the estimator variance from all the components. 

The signal is expected to scale as Rf and the noise as R^ 2 in a situation where 
the baseline distribution is uniform ie. Pn(U,v) is independent of U. This holds at 
U < 300 for the GMRT (Figure I3~5l) . and the expected scaling is seen for Rf, > 20 Mpc. 
For smaller bubbles the signal extends to larger baselines where pn{U,v) falls sharply, 
and the signal and the noise both have a steeper R b dependence. The MWA baseline 
distribution is flat for only a small U range (Figure 13. 5[) beyond which it drops. In this 
case the signal and noise are found to scale as Rf and R\ respectively. Note that the 
maximum baseline at MWA is U — 750, and hence a considerable amount of the signal 
is lost for R b < 10 Mpc. 

At both the GMRT and the MWA, for 100 hrs of observation, the noise is larger than 
the signal for bubble size Rb < 40 Mpc. At the other extreme, for an integration time 
of 10 5 hrs the noise is below the signal for R b > 6 Mpc for the GMRT and R b > 8 Mpc 
for the MWA. The foreground contribution turns out to be smaller than the signal for 
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Figure 3.7: Same as Figure I3T61 except that Filter II is used instead of Filter I. 
the entire range of bubble sizes that we have considered, thus justifying our choice of 
filters. Note that Filter II is more efficient in foreground subtraction, but it requires prior 
knowledge about the frequency dependence. For both the filters the foreground removal 
is more effective at the GMRT than the MWA because of the frequency dependence of 
Pn{U, v). The assumption that this is proportional to A 2 is valid only when p(U, v) is 
independent of U, which, as we have discussed, is true for a large U range at the GMRT. 
The A dependence is much more complicated at the MWA, but we have not considered 
such details here as the foreground contribution is anyway smaller than the signal. It 
should also be noted that the foreground contribution increases at small baselines (eq. 
13.101) . and is very sensitive to the smallest value of U which we set at U = 20 for our 
calculations. Here it must be noted that our results are valid only under the assumption 
that the foregrounds have a smooth frequency dependence. A slight deviation from this 
and the signal will be swamped by the foregrounds. Also note that this filtering method 
is effective only for the detection of the bubbles and not for the statistical HI fluctuations 
signal. 

The contribution from the HI fluctuations impose a lower limit on the size of the 
bubble which can be detected. However long be the observing time, it will not be 
possible to detect bubbles of size Rb < 8 Mpc using the GMRT and size Rb < 16 Mpc 
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using the MWA. The HI fluctuation contribution increases at small baselines. The 
problem is particularly severe at MWA because of the dense sampling of the small 
baselines and the very large field of view. We note that the MWA is being designed 
with the detection of the statistical HI fluctuation signal in mind, and hence it is not 
surprising that this contribution is quite large. For both telescopes it may be possible to 
reduce this component by cutting off the filter at small baselines. We have not explored 
this possibility in this work because the enormous observing times required to detect 
such small bubbles makes it unfeasible with the GMRT or MWA. 

Figure 13.81 shows the observation time that would be required to detect bubbles of 
different sizes using Filter I for GMRT and the MWA. Note that the observing time 
shown here refers to a 3cr detection which is possibly adequate for targeted searches 
centered on observed quasar position. A more stringent detection criteria at the 5 a 
level would be apropriate for a blind search. The observing time would go up by a factor 
of 3 for a 5 a detection. The observing time is similar for Filter II and hence we do not 
show this separately. In calculating the observing time we have only taken into account 
the noise contribution as the other contributions do not change with time. The value of 
Rb below which a detection is not possible due to the HI fluctuations is shown by vertical 
lines for both telescopes. We see that with 100 hrs of observation both the telescopes 
will be able to detect bubbles with > 40 Mpc while bubbles with Rt> > 22Mpc can 
be detected with 1000 hrs of observation. 

The possibility of detecting a bubble is less when the bubble centre does not coincide 
with the centre of the field of view. In fact, the SNR falls as e~ ^ e ° if the bubble center 
is shifted away by 6 C from the center of the field of view and the filter is also shifted so 
that its center coincides with that of the bubble. There will be a corresponding increase 
t b s oc e 29 ^ o in the observing time required to detect the bubble. It will be possible to 
detect bubbles only if they are located near the center of the field of view (& c <C Q ), and 
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Figure 3.8: The observing time t Q b s that would be required for a 3 cr detection of a bubble 
of radius R b provided it is at the center of the field of view. The vertical 
lines shows the lower limit (due to HI fluctuations) where a detection will be 
possible (R b = 8 Mpc for GMRT and R b = 16 Mpc for MWA). 
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Figure 3.9: The Overlap between the signal and the filter when there is a mismatch AO 
between the centers of the bubble and the filter for GMRT (left) and MWA 
(right). The results are shown for different bubble sizes. 

the required observing time increases rapidly with 8 C for off-centered bubbles. 

When searching for bubbles in a particular observation it will be necessary to consider 
filters corresponding to all possible value of Rb, v c and 6 C . A possible strategy would 
be to search at a discrete set of values in the range of Rb, v c and 8 C values where a 
detection is feasible. The crucial issue here would be the choice of the sampling density 
so that we do not miss out an ionized bubble whose parameters do not exactly coincide 
with any of the values in the discrete set and lie somewhere in between. To illustrate 
this we discuss the considerations for choosing and optimal value of A8 C the sampling 
interval for 9 C . We use (E)[A9] to denote the expectation value of the estimator when 
there is a mismatch A9 between the centers of the bubble and the filter. The ratio 
Overlap = (E)[A0]/(E)[0], shown in Figure EH for GMRT (left panel) and MWA (right 
panel), quantifies the overlap between the signal and the filter as A9 is varied. We see 
that the choice of A8 would depend on the size of the bubble we are trying to detect and 
it would be smaller for the GMRT as compared to the MWA. Permitting the Overlap 
to drop to 0.9 at the middle of the sampling interval, we find that it is 8' at the GMRT 
and 20' at the MWA for R b = 50 Mpc. 

The MWA is yet to be constructed, and it may be possible that an antennae distribu- 
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Figure 3.10: Same as the Figure 13.81 considering three different antenna distributions 
Pantir) oc 1/r, 1/r 2 , 1/r 3 for the MWA. 
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tion different from p an t{ r ) 1/r 2 , may improve the prospects of detecting HII bubbles. 
We have tried out p an t{f) 1/r and 1/r 3 for which the results are shown in Figure 13.101 
We find that the required integration time falls considerably for the 1/r 3 distribution 
whereas the opposite occurs for 1/r. For example, for Rb = 50 Mpc the integration time 
increases by 5 times for 1/r and decreases by 3 times for 1/r 3 as compared to 1/r 2 . 
Based on this we expect the integration time to come down if the antenna distribution 
is made steeper, but this occurs at the expense of increasing the HI fluctuations and the 
foregrounds. We note that for the 1/r 3 distribution the foreground contribution is more 
than the signal, but it may be possible to overcome this by modifying the filter. The 
increase in the HI fluctuations is inevitable, and it restricts the smallest bubble that 
can be detected to R = 26MPc for 1/r 3 . In summary, the 1/r 2 distribution appears 
to be a good compromise between reducing the integration time and increasing the HI 
fluctuations and foregrounds. 

Finally we examine some of the assumptions made in this work. First, the Fourier 
relation between the specific intensity and the visibilities (eq. 13. 1 p will be valid only near 
the center of the field of view and full three dimensional wide-field imaging is needed 
away from the center. As the feasibility of detecting a bubble away from the center falls 
rapidly, we do not expect the wide-field effects to be very important. Further, these 
effects are most significant at large baselines whereas most of the signal from ionized 
bubbles is in the small baselines. 

Inhomogeneities in the IGM will affect the propagation of ionization fronts, and 



the i onized bubbles are not expected to be exactly spherical (jWyithe. Loeb fc Barnes 



20051 ). This will cause a mismatch between the signal and the filter which in turn will 
degrade the SNR. In addition to this, in future we plan to address a variety of other 
issues like considering different observing frequencies and making predictions for the 
other upcoming telescopes. 
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Terrestrial signals from television, FM radio, satellites, mobile communication etc., 
collectively referred to as RFI, fall in the same frequency band as the redshifted 21cm 
signal from the reionization epoch. These are expected to be much stronger than the 
expected 21cm signa 



Bowman et al. 



and it is necessary to quantify and characterize the RFI. Recently 
(120071 ) have characterized the RFI for the MWA site on the frequency 
range 80 to 300MHz. They find an excellent RFI environment except for a few channels 
which are dominated by satellite communication signal. The impact of RFI on detecting 
ionized bubbles is an important issue which we plan to address in future. 

The effect of polarization leakage is another issue we postpone for future work. This 
could cause polarization structures on the sky to appear as frequency dependent ripples 
in the foregrounds intensity . This could be particularly severe for the MWA. 
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4 Simulating Matched Filter Search 
for Ionized Bubbles^ 



4.1 Introduction 

In the previous Chapter we present an analytic framework for predicting the expected 
value and the standard deviation o of the matched filter estimator for the detection of a 
spherical ionized bubble of comoving radius Rf,. We identify three different contributions 
to a, namely foregrounds, system noise and the fluctuations in the HI outside the bubble 
that we are trying to detect. Our analysis shows that the matched filter effectively 
removes the foreground contribution so that it falls below the signal. Considering the 
system noise for the GMRT and the MWA we find that a 3 a detection will be possible for 



This Chapter is adapted from the paper "Simulating the impact of HI flu ctuations on 



matched filter search for ionized bubbles in redshifted 21 cm maps" by \Datta et al. 
\200A ). 
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a bubble of comoving radius R b > 40 Mpc in 100 hrs of observation and R b > 22 Mpc in 
1000 hrs of observation for both the instruments. The HI fluctuations, we find, impose a 
fundamental restriction on bubble detection. Under the assumption that the HI outside 
the ionized bubble traces the dark matter we find that it is not possible to detect bubble 
of size R b < 8 Mpc and R b < 16 Mpc at the GMRT and MWA respectively. Note that 
the matched filter technique is valid for both, a targeted search around QSOs as well as 
for a blind search in a random direction. 

Here we validate the visibility based matched filter technique introduced in the pre- 
vious Chapter through simulations of bubble detection. Our simulations are capable 
of handling interferometric arrays with widely different configurations like the GMRT 
and the MWA , the two instruments that we consider here. As mentioned earlier, the 
fluctuations in the HI outside the target bubble impose a fundamental restriction for 
bubble detection. The analytic approach of the Chapter 3 assumes that the HI outside 
the bubble traces the dark matter. In this Chapter we carry out simulations that incor- 
porate this assumption and use these to assess the impact of HI fluctuations for bubble 
detection. We also use the simulations to determine the accuracy to which the GMRT 
and the MWA will be able to determine the size and the position of an ionized bubble, 
and test if this is limited due to the presence of HI fluctuations. In a real situation a 
typical FoV is expected to contain several ionized patches besides the one that we are 
trying to detect. We use simulations to assess the impact of HI fluctuations for bubble 
detection in patchy reionization scenarios. 

The outline of the Chapter is as follows. Section 4.2 presents a brief description of how 
we simulate 21-cm maps for three different scenarios of the HI distribution, one where the 
HI traces the dark matter and two with patchy reionization. Subsections 4.2.1 and 4.2.2 
respectively discuss how the simulated maps are converted into visibilities and how the 
matched filter analysis is simulated. We present our results in Section 4.3. Subsections 
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4.3.1, 4.3.2 and 4.3.3 present results for bubble detectability, size determination and 
position determination under the assumption that the HI outside the bubble traces the 
dark matter. Section 4.3.4 presents results for bubble detectability in patchy reionization 
scenarios. We discuss redshift dependence of bubble detection in Section 4.4 and present 
our summary in Section 4.5. 



For the GMRT we have used the telescope 
the MWA we use the telescope parameters from 



parameters from their website, while for 



Bowman et al. 



(120061 ) . The cosmological 



parameters for a flat (k = 0) ACDM model used throughout this paper are Q r , 
0.3, tt b h 2 = 0.022, n s = l.,h = 0.74, <r 8 = 1. 



4.2 Method of Simulation 



We have simulated the detection of the HI signal of an ionized bubble whose center is 
at redshift z c = 6 which corresponds to v c = 203 MHz. The choice of z value is guided 
by the fact that we expect large ionized regions towards t h e end of reionization z > 6 



(IWvithe fc Loeb 



20041 : 



Furlanetto. McQuinn &: Hernquist 



20061 ). Our aim here is to 



validate the analytic calculations of the Chapter 3 and hence the exact value of z is not 
very important. 

We consider four scenarios of reionization for bubble detection. In the first three 
scenarios there is a spherical ionized bubble, the one that we are trying to detect, at 
the center of the FoV. This bubble has comoving radius Rb and is embedded in HI that 
traces the dark matter. In the first scenario there is a single bubble in the field of view. 
We refer to this as the SB scenario. In this scenario the HI fraction xm is assumed to be 
uniform outside the bubble. The uncertainty due to the HI fluctuations is expected to be 
lowest in this scenario because of the absence of patchiness. This is the most optimistic 
scenario for bubble detection. 
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In the next two scenarios, we attempt to quantify the effect of patchy reionization 
(PR) outside the bubble that we are trying to detect by introducing many other, pos- 
sibly overlapping, bubbles in the FoV. Unfortunately, there is no obvious way to fix 
the sizes of these bubbles from any theoretical models as they depend crucially on the 
nature of reionization sources and other physical factors. In scenario PR1, we assume 
that the large HII regions which we are trying to detect are surrounded by many small 
ionized regions whose sizes are fixed by the following procedure: we assume the globally 
averaged neutral fraction xhi to be ~ 0.5; the reason for this choice is that the effects of 
patchiness would be most prominent when typically half of the IGM is ionized. Given 
the value of xhi, we try to obtain a reasonable estimate of the size of the background 
rom available models. For ex ample, semi- numeric simulations of patchy reion- 



bubbles 

ization (jMesinger fc Furlanettol . 120071 ) predict that the bubble size distribution peaks 
around 5Mpc when xhi = 0.61 (see their Fig 6). We thus choose the spherical back- 
ground bubbles to have radii 6 Mpc and compute the number of background bubbles by 
demanding that the resulting neutral fraction is 0.5. The bubble centres are chosen such 
that they trace the underlying dark matter distribution. At the end, the value of xm 
turns out to be slightly higher 0.62 because of overlap of the bubbles. Note that because 
of these overlaps, the shapes of the resulting ionized regions would not always be per- 
fectly spherical. In this scenario, we have essentially attempted to capture a situation 
where there are many small, possibly overlapping ionized regions produced by galaxies 
and a few large ionized regions (like the one that we are trying to detect) produced by 
QSOs. 

Since the choice of the background bubble size is not robust by any means, we 
consider a different scenario PR2 where these bubbles have the same comoving radius 
Rb as the bubble that we are trying to detect. The centers of these extra bubbles trace 
the dark matter distribution as in PR1. The number of bubbles is fixed by the globally 
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averaged xhi which we take to be 0.62 same as in PR1. The PR2 scenario represents 
a situation where we predominantly have large ionized regions produced either by rare 
luminous sources or through the overlap of several small ionized regions in the later 
stages of reionization. 

A particle-mesh (PM) N-body code was used to simulate the dark matter distribution. 
In the Chapter 3 we show that the HI signal of the ionized bubble is largely concentrated 
at small baselines or large angular scales, thus a very high spatial resolution is not 
required. We have used a grid spacing of 2 Mpc for the simulations. This is adequate 
for bubbles in the range 4 < < 50 Mpc that we consider. The simulations use 256 3 
particles on a 256 3 mesh. For the GMRT a single N-body simulation was cut into 8 equal 
cubes of size 256 Mpc on each side. Considering that each cube may be viewed along 
three different directions, we have a total of 24 different realization of the dark matter 
distribution. Each cube corresponds to 18 MHz in frequency and ~ 2° in angle which 
is comparable to the GMRT FoV which has FWHM=1.7° at 203MHz. The MWA FoV 
is much larger (FWHM=13°). Here eight independent N-body simulations were used. 
Viewing these along three different directions gives twenty four different realizations of 
the dark matter distribution. Limited computer memory restricts the simulation size 
and the angular extent (~ 4°) is considerably smaller than the MWA FoV. We do not 
expect this to affect the signal but the contribution from the HI fluctuations outside the 
bubble is possibly underestimated for the MWA. 

The dark matter density contrast 5 was used to calculate the redshifted 21 — cm 
specific intensity I v = I u xm(l + 8) for each grid point of our simulation. Here I u = 
2.5 x 10 2 ^ (^P) (x) (iff})) and Xm tne hydrogen neutral fraction is inside the 
ionized bubbles and 1 outside. The simulated boxes are transformed to frequency and 
sky coordinate. Figure 14.11 shows the HI image on a slice through the center of the 
bubble of radius Rb = 20 Mpc. The mean neutral fraction xhi is ~ 1 in the SB scenario, 
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Figure 4.1: This shows HI images on slices through the center of the bubble for the 
four different scenarios SB, PR1, PR2 and SM (from left to right). In first 
three panels the central, circular dark region of radius Rb = 20 Mpc shows 
the HII bubble that we are trying to detect. The HI outside this bubble 
traces the dark matter distribution. In the SB scenario (left) the hydrogen 
neutral fraction is a; hi = 1 outside the bubble. In the PR1 scenario (2nd 
from left) the extra bubbles are all of a fixed comoving radius 6 Mpc. In 
the PR2 scenario (3rd from left) the extra bubbles have the same comoving 
radius as the bubble that we are trying to detect. In both the PR1 and PR2 
reionization scenarios the centers of the extra bubbles trace the dark matter 
distribution and % = 0.62. In the SM scenario (right) the central region 
up to radius 27 Mpc is fully ionized (marked with solid circle) and beyond 
that region up to radius 42 Mpc the region is partially filled with HI patches 
(dashed circle). The mean neutral fraction is xm = 0.5. These simulations 
are all for the GMRT. 

while it is ~ 0.62 for the two PR simulations shown here. 

The three scenarios discussed above consider only spherical bubbles, and the only 
departures from sphericity arise from bubble overlap. It is important to assess how 
well our bubble detection technique works for non-spherical bubbles, which we do using 
ionization maps produced by the semi-numeric (SM) approach. In part icular, we use 



maps obtained by the method of 



Choudhury. Haehnelt fc Regan 



(120081 ). Essentially, 



these maps are produced by incorporating an excursion-set based technique for identify- 



ing ionized regions given the 



2007 



Mesinger &: Furl anetto. 



density distribution and the ionizing sources (IZahn et al. 



2007 



Choudhury. Haehnelt fc Regan 



Geil fc Wvithe 



20081 ). In addition, the method of 



(120081 ) incorporate inhomogeneous recombination and 



self-shielding of high-density regions so that it is consistent with the "photons-starved" 
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reionization scenario implied by the Lyo forest data (IBolton &: Haehnelt 



2007 



Choudhury. Ferrara fc 



20081 ). We use a simulation box of size 270 Mpc with 2000 3 particles which can resolve 
collapsed halos as small as ~ 10 9 M Q . The ionization maps are generated at a much lower 
resolution with a grid size of 2.7 Mpc. The box corresponds to 19 MHz in frequency 
and ~ 2° in angle comparable to the GMRT FoV. We have assigned luminosities to the 
collapsed halos such that the mean neutral fraction xm = 0.5. The most massive halo 
(mass ~ 1O 13 M ) identified in the box is made to coincide with the box centre and we 
assume that it hosts a luminous QSO; its luminosity and age are chosen such that it 
would produce a spherical HII region of comoving si ze ~ 27 Mpc in a com pletely ho- 



Geil fc Wvithel ((20081)]. However 



mogeneous neutral medium [see, e.g., equation (8) of 
the actual ionized region is far from spherical both because of the surrounding bubbles 
from other halos and also because of inhomogeneous recombination. We find visually 
from the maps (see the rightmost panel of Figure fl~Tl) that the HII region is fully ionized 
up to radius ~ 27 Mpc. Beyond that the region is partially filled with neutral patches. 
This patchy ionized region extends up to radius ~ 42 Mpc and then merges with the 
average IGM. The fully ionized region and the region with HI patches are marked with 
two circles. We use this box for GMRT as three independent realizations viewing the 
box along three different directions. For the MWA we need a much larger simulation 
box which requires substantially more computing power, beyond the resources available 
to us at present. Hence we do not consider the MWA for this scenario. 
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4.2.1 Simulating visibilities 



The quantity measured in radio-interferometric observations is the visibility V(U, v) 
which is related to the specific intensity pattern on the sky I u (0) as 

V(U, u) = J d 2 9A(9)I u (9)e 2mlV (4.1) 

Here the baseline U = d/X denotes the antenna separation d projected in the plane 
perpendicular to the line of sight in units of the observing wavelength A, 9 is a two 
dimensional vector in the plane of the sky with origin at the center of the FoV, and 
A(9) is the beam pattern of the individual antenna. For the GMRT this can be well 
approximated by Gaussian A{9) = e~ e2 ^ e ° 2 where #0 « 0.6 #fwhm an d we use the values 
1.7° for # at 203 MHz corresponding to the redshift z = 6 for the GMRT. The MWA 
beam pattern is expec ted to be quite comp licated, and depends on the pointing angle 



relative to the zenith (IBowman et al. 



20071 ). Our analysis largely deals with the beam 



pattern within 2° of the pointing angle where it is reasonable to app roximate the beam 



Bowman et al 



20071 ) . We approximate 



as being circularly symmetric (Figures 3 and 5 of [ 
the MWA antenna beam pattern as a Gaussian. 

We consider 128 frequency channels across 18MHz bandwidth. The image I v {9) at 
each channel is multiplied with the telescope beam pattern A(8, v). The discrete Fourier 
transform (DFT) of the product I u (9) A(9, v) gives the complex visibilities V(U, v). The 
GMRT simulations have baselines in the range 30.5 < U < 3900 which is adequate to 
capture the HI signal from ionized bubbles which is expected to be confined to small 
baselines U < 1000. 

The visibility recorded in radio-interferometric observations is actually a combination 
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Figure 4.2: This shows the visibility signal (real part) from a frequency slice through the 
center of a spherical ionized bubble of comoving radius 20 Mpc embedded in 
HI. The solid curves show the expected signal assuming that the bubble is 
embedded in uniformly distributed HI. The data points show the visibilities 
for a few randomly chosen baselines from our simulation of the SB scenario. 
The difference between the data points and the solid curve is due to the 
fluctuations in the HI outside the bubble. Each panel corresponds to a 
different realization of the simulation. 
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Figure 4.3: Same as the previous figure except that U is fixed at 110 while the frequency 
varies , and Au = v — u c . 

of several contributions 

V(U, v) = S(U, v) + HF(U, v) + N(U, v) + F(U, v) . (4.2) 

where S(U, v) is the HI signal that we are interested in, HF(U, v) is contribution from 
the fluctuating HI outside the bubble that we are trying to detect, N(U, v) is the system 
noise which is inherent to the measurement and F(U, v) is the contribution from other 
astrophysical sources referred to as the foregrounds. The signal S(U, v) from an ionized 
bubble of comoving radius Rb embedded in an uniform HI distribution can be analytically 
calculated (Chapter 3). The solid curve in Figures 14721 and I4~3l show the expected signal 
for Rb = 20Mpc. The U extent, frequency extent and peak value of the signal scale as 
Rk , Rb and Rl respectively for other values of Rb- Note that S(U, v) is real when the 
bubble is at the center of the FoV. 

The data points shown in Figures 14.21 and 14.31 are the real part of a few randomly 
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chosen visibilities determined from the simulation of a Rj, = 20 Mpc bubble in the SB 
scenario. The deviations from the analytic predictions are due to the HI fluctuations 
HF(U, v) ie. in the SB scenario the HI outside the bubble traces the dark matter 
fluctuations. Notice that these fluctuations are often so prominent that the signal cannot 
be made out. We expect even larger fluctuations in the other three scenarios which 
incorporate patchiness of reionization. 

The system noise contribution N(U, v) in each baseline and frequency channel is 
expected to be an independent Gaussian random variable with zero mean ((N) = 0) 
and variance \J (iV 2 )is independent of U and v c . We use (Chapter 3) 

where C x has values 0.53Jy and 54.21Jy for the GMRT and the MWA respectively 
(Chapter 3). 

The contribution from astrophysical foregrounds F(U, v) is expected to be several 
order of magnitude stronger than the HI signal. The foregrounds are predicted to have a 
featureless, continuum spectra whereas the signal is expected to have a dip at u c (Figure 
14. 3D . This difference holds the promise of allowing us to separate the signal from the 
foregrounds. 



4.2.2 Simulating signal detection 

The signal component S(U, v) in the observed visibilities V(U, v) is expected to be buried 
deep in other contributions many of which are orders of magnitude larger. Detecting 
this is a big challenge. For optimal signal detection we consider the estimator (Chapter 
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3) 

E = Y, S f(U a ,Vb)V(Ua,Vb) (4-4) 

a,b 

where Sf(U, v) is a filter which has been constructed to detect a particular ionized 
bubble, V{U a ,Vb) refer to the observed visibilities and U a and refer to the different 
baselines and frequency channels in the observation. The filter Sf(U,v) depends on 
[Rf, z c , 9 C ] the comoving radius, redshift and angular position of the bubble that we are 
trying to detect. We do not show this explicitly, the values of these parameters will be 
clear from the context. 

The baselines obtained using DFT in our simulations are uniformly distributed 
on a plane. In real observations, the baselines will have a complicated distribution 
depending on the antenna layout and direction of observation. We incorporate this 
through the normalized baseline distribution function p^(U,v) which is defined such 
that cPUdv pn(U,v) is the fraction of data points ie. baselines in the interval d 2 U dv 
and J d 2 U J dvp N (U,is) = 1. We use the functional forms of pn determined in the 
Chapter 3 for the GMRT and the MWA. 

Using the simulated visibilities, we evaluate the estimator as 

E = (AU) 2 Av S *f(Ua, v b )V(Ua, Vb)p N (U a , (4.5) 

a,b 

where the sum is now over the baselines and frequency channels in the simulation. 
The filter Sf(U, v) (Filter I of Chapter 3) is defined as 

\ 2 

v 



S f (U,u)=l-\ S(U,u) 



e ( i-2|„-„j/ B *) r +B 'v y)dl/ 



B 



v c -B'/2 
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(4.6) 

where the first term S(U, v) is the expected signal of the bubble that we are trying 
to detect. We note that this term is the matched filter that gives the maximum signal 
to noise ratio (SNR). The second term involving the Heaviside function 0(x) subtracts 
out any frequency independent component from the frequency range v c — B' /2 to u c + 
B' /2. The latter term is introduced to subtract out the foreground contributions. The 
(z//z/ c ) 2 term accounts for the fact that p N (U,v) changes with frequency (equivalently 
wavelength) . 

We have used the 24 independent realizations of the simulation for the first three 
scenarios to determine the mean (E) and the variance ((AE) 2 ) of the estimator. The 
high computational requirement restricts us to use just 3 realizations for the SM scenario. 
Only the signal is correlated with the filter, and only this is expected to contribute to 
the mean (E). All the other components are uncorrelated with the filter and they are 
expected to contribute only to the variance ((AE) 2 ). The variance is a sum of three 
contributions (Chapter 3) 

((AE) 2 ) = ({AE) 2 ) HF + ((AE) 2 ) N + ((AE) 2 ) pc . (4.7) 

The simulations give an estimate of ((AE) 2 )hf the contribution from HI fluctuations. 
We do not include system noise explicitly in our simulations. The noise contribution 
from a single visibility (eq. 14.31) is used to estimate ((AE) 2 ) N (eq. 3.19 of Chapter 
3). Under the assumed foreground model, the foreground contribution ((AE) 2 ) FG is 
predicted to be smaller than the signal and we do not consider it here. 
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4.3 Results 

We first consider the detection of an ionized bubbles and the estimation of its parameters 
in the SB scenario where there is only a single bubble in the FoV. We consider the most 
optimistic situation where the bubble is located in the center. In reality this can only be 
achieved in targeted observations of ionized bubbles around luminous QSOs. In a blind 
search, the bubble in general will be located at some arbitrary position in the FoV, and 
not the center. It has already been mentioned that the foregrounds can be removed by 
a suitable choice of the filter. Further, the system noise can, in principle, be reduced 
by increasing the observation time. The HI fluctuations outside the bubble impose a 
fundamental restriction on bubble detection. 

4.3.1 Restriction on bubble detection 

We have carried out simulations for different values of the bubble radius Rb chosen 
uniformly at an interval of 2 Mpc in the range 4 to 50 Mpc. In each case we consider 
only the most optimistic situation where the bubble radius Rf used in the filter is 
precisely matched to Rb- In reality it is necessary to try filters of different radius Rf to 
determine which gives the best match. 

Figures 14.41 and 14.51 shows the results for the GMRT and the MWA respectively. We 
compare the analytic predictions of the Chapter 3 (left panel) with the prediction of 
our simulations (right panel). The analytical predictions for the mean value (E) arising 
from the signal and \J ((AE) 2 ) ^ due to the HI fluctuations are respectively calculated 
using equations (3.15) and (3.22) of the Chapter 3. The signal depends on the bubble 
radius Rb and the mean neutral fraction which is taken to be xhi = 1. The uncertainty 
due to the HI fluctuations is calculated using the dark matter power spectrum under 
the assumption that the HI traces the dark matter. 
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Figure 4.4: The estimator E (defined in equation I4.4D for bubble size Rb ranging from 
4 Mpc to 50 Mpc for the GMRT in the SB scenario. It is assumed that the 
filter is exactly matched to the bubble. The left panel shows the analytic 
predictions for the mean estimator (E) and the 3 — a error-bars due to the HI 
fluctuations. The solid and the dashed lines in the right panel respectively 
show the (E) and the 3 — a envelope determined from the simulations. The 
data points in the right panel show E in the individual realizations. 
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We find that (E) and y ((AE) 2 } determined from the simulations is in rough 
agreement with the analytic predictions. The mean (E) is in very good agreement for 
R b > 6Mpc, there is a slight discrepancy for smaller bubbles arising from the finite 



grid size (2Mpc in the simulation). The HI fluctuations y ((AE) 2 )^^ are somewhat 
underestimated by the simulations. This is more pronounced for the MWA where the 
limited box size of our simulations results in a FoV which is considerably smaller than 
the actual antennas. We note that the 24 different values of E determined from the 



different realizations of the simulation all lie within (E) ± 3y {(AE) 2 ) determined 
from the analytic predictions. 

The good agreement between the simulation results and the analytical predictions is 
particularly important because each is based on several approximations, many of which 
differ between the two methods. Our results show that the effect of these approximations, 
though present, are well under control. The analytical method has the advantage that 
it is very easy to calculate and can be evaluated very quickly at an extremely low 
computational cost. Unfortunately, its utility is mainly limited to the SB scenario and it 
cannot be easily applied to an arbitrary PR scenario with a complicated HI distribution. 
Simulations, though computationally more cumbersome and expensive, are useful in 
such a situation. It is thus important to test that the two methods agree for the SB 
scenario where both of them can be applied. Note that the HI fluctuation predicted by 
the SB scenario sets the lower limit for the HI fluctuation in any of the PR models. It 
is expected that patchiness will increase the HI fluctuations above the SB predictions. 

It is meaningful to attempt bubble detection at, say 3a confidence level, only if 
(E) > 3\J ((AE) 2 )hf- The HI fluctuations overwhelm the signal in a situation where 
this condition is not satisfied, and bubble detection is not possible. In a situation where 
this condition is satisfied, an observed value E of the estimator can be interpreted as 
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a 3<j detection if E Q > 3y ((AE) 2 }. The simulations show that a 3 — a detection is not 
possible for R b < 6 Mpc and R b < 12 Mpc at the GMRT and MWA respectively. As 
noted earlier, the HI fluctuations are somewhat under predicted in the simulations and 
the analytic predictions Rb < 8 Mpc and Rb < 16 Mpc respectively, are somewhat larger. 

The limitation on the bubble size Rb that can be detected is larger for the MWA as 
compared to the GMRT. This is because of two reasons, the first being the fact that 
the MWA has a very dense sampling of the small baselines where the HI fluctuation 
are very large, and the second being the large FoV. In fact, the baseline distribution of 
the experiment has a significant role in determining the quantum of HI fluctuations and 
thereby determining the lower cut-off for bubble detection. Looking for an optimum 
baseline distribution for bubble detection is also an issue which we plan to address in 
future. In a situation where the antenna layout is already in place, it may possible to 
tune the filter to reduce the HI fluctuations. 



We have not considered the effect of peculiar velocities (jBharadwaj fc Ali . 



20M) in 



our simulations. From equation (3.22) of the Chapter 3 we see that the HI fluctua- 



tions scale as y ((AE) 2 } HF oc where C\ is the HI multi-frequency angular power 
spectrum (MAPS). The C/s increase by a factor ~ 2 due to peculiar velocities, whereby 



((AE) 2 ) hf goes up by a factor ~ 1.5. This increase does not significantly change our 
results, and is small compared to the other uncertainties in the PR models. 



The signal (E) and the HI fluctuations y ((AE) 2 )hf both scale as oc xhi, and the 
lower limit for bubble detection is unchanged for smaller neutral fractions. 



4.3.2 Size determination 

In this subsection we estimate the accuracy to which it will be possible to determine the 
bubble radius Rb- This, in general, is an unknown quantity that has to be determined 
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from the observation by trying out filters with different values of Rf. In the matched 
filter technique we expect the predicted SNR (only system noise) ratio 

IE) 

SNR = V I (4.8) 
\A(A£) 2 >ns 

to peak when the filter is exactly matched to the signal ie. Rf — Rb- The solid line 
in the right panel of Figure 14.61 shows this for Rf, = 10 Mpc. We find that the SNR 
peaks exactly when the filter size Rf = 10 Mpc. We propose that this can be used to 
observationally determine Rj,. For varying Rf, we consider the ratio of the observed 
value E Q to the expected system noise y ((A£') 2 )ns, referring to this as the SNR. The 
Rf value where this SNR peaks gives an estimate of the actual bubble size Rb- The 
observed SNR will differ from the predictions of eq. (14. 8p due to the HI fluctuations 
outside the bubble. These variations will differ from realization to realization and this 
can introduce uncertainties in size estimation. We have used the simulations to estimate 
this. 

The left panel of Figure l4~o! shows the SNR as a function of Rf for 4 different real- 
izations of the simulation for the GMRT with bubble size Rb = 10 Mpc. We see that for 
Rf < Rb the SNR shows a very similar behavior in all the realizations, and it always 
peaks around 10 Mpc as expected. For Rf > R^ the behavior of the SNR as a function of 
Rf shows considerable variation across the realizations. In some cases the drop in SNR 
away from the peak is quite rapid whereas in others it is very gradual (for example, the 
dashed-dot-dot curve). In many cases there is an spurious extra peak in the SNR at an 
Rf value that is much larger than Rb. These spurious peaks do not pose a problem for 
size determination as they are well separated from R^ and can be easily distinguished 
from the genuine peak. 

The error-bars in the right panel of Figure 14.61 show the 3 — a fluctuation in the 
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Figure 4.6: The SNR = (E) / J ((AE) 2 ) NS for 1000 hrs observation with the GMRT as 
a function of the filter size Rf for the case where the actual bubble size is 
Rb = 10 Mpc . The left panel shows 4 different realizations of the simulation. 
The right panel shows the mean SNR and 3 — a error-bars calculated using 
24 realizations. The solid line shows the analytical predictions. 
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simulated SNR determined from 24 realizations of the simulation. Note that the fluc- 
tuations at different Rf are correlated. Although the overall amplitude changes from 
one realization to another, the shape of the curve in the vicinity of Rf = Rb is nearly 
invariant across all the realizations. In all of the 24 realizations we can identify a well 
defined peak at the expected value Rf = Rb. 

Figures 14.71 and 14.81 show the results for a similar analysis with Rb = 20 Mpc for 
the GMRT and the MWA respectively. It is not possible to detect a bubble of size 
Rb = 10 Mpc with the MWA, and hence we do not show this. Here again,we find that 
for all the realizations of the simulations the SNR peaks at Rf = Rb- The relative 
variations in the SNR across the realizations is much less for Rb = 20 Mpc as compared 
to 10 Mpc and there are no spurious peaks. Also, for the same bubble size the variations 
are smaller for the GMRT as compared to the MWA. We do not find any spurious peaks 
for R b = 20 Mpc. 

A point to note is that the mean SNR determined from the simulations is somewhat 
smaller than the analytic predictions, both being shown in the right panels of Figures 
14.61 14. 71 and Figure I4~%1 There are a couple of reasons that could account for this namely, 
(i) the bubble in the simulation is not exactly a sphere because of the finite grid size and 
thus the match between the filter and the signal is not perfect even when the sizes are 
same and (ii) the finite box-size imposes a minimum baseline beyond which the signal 
is not represented in the simulation. 

Based on our results we conclude that in the SB scenario for the GMRT the accuracy 
to which the bubble size can be determined in our simulations is decided by the resolution 
2 Mpc and not by the HI fluctuations. In reality the limitation will come from the angular 
resolution of the instrument which sets the limit at 0.5 Mpc for the GMRT and 8 Mpc 
for the MWA. 
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Figure 4.7: Same as the Figure WE for R b = 20 Mpc for the GMRT. 

The height of the SNR peak depends on the neutral fraction and it can be used 

to observationally determine this. We find that the HI fluctuations do not change the 

position of the peak but introduce considerable variations in its height even if xm — 1- 

The HI fluctuations restrict the accuracy to which the neutral fraction can be estimated, 

an issue that we propose to address in future work. 



4.3.3 Determining the position 

In the previous two subsections, we have considered cases where the bubble's position is 
known. Here we assume that the bubble's size is known and we estimate the accuracy to 
which its position can be determined in the presence of HI fluctuations. The situation 
considered here is a blind search whereas the former is a targeted search centered on a 
QSO. 

In a real situations it would be necessary to jointly determine four parameters the 
bubble radius Rb, two angular coordinates (9 x ,9 y ) and the central frequency v c from 
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Figure 4.8: Same as the Figure HSl for R b = 20 Mpc for the MWA. 
the observation. However, to keep the computational requirement under control, in 
this analysis we assume that Rb is known. The bubble is placed at the center of the 
FoV and frequency band, and we estimate how well the position can be recovered from 
the simulation. To determine the bubble's position we move the center of the filter 
to different positions and search for a peak in the SNR. To reduce the computational 
requirement, this is done along one direction at a time, keeping the other two directions 
fixed at the bubble's actual center. We have also carried out simulations where the 
bubble is located off-center. We do not explicitly show these results because they are 
exactly the same as when the bubble is at the center except for the fact that the value 
of the peak SNR is lower because of the primary beam pattern. 

Figure I4.9I shows the results for Rb = 10 Mpc for the GMRT. The left panel shows 
results for 4 realizations of the simulation, the right panels show the mean and 3 — a 
determined from 24 realizations of the simulations and the analytic prediction for the 
mean value. In all peak is seen at the expected position matched with the 

bubble's actual center. The HI fluctuations pose a severe problem for determining the 
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Figure 4.9: The SNR for 1000 hrs GMRT observations for a bubble of size R b = lOMpc 
located at the center of the filed of view. The filter scans along 6 X , 8y, Vf 
(top, middle, bottom) to determine the bubble's position. The left panel 
shows results for 4 realizations of the SB simulation, the right panels show 
the mean (dashed curve) and 3 — a error-bars determined from 24 realizations 
of the simulation and the analytic prediction for the mean (solid curve). 
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bubble's position as it introduces considerable fluctuations in the SNR. In some cases 
these fluctuations are comparable to the peak at the bubble's actual position (see the 
dashed line in the upper left panel). The possibility of spurious peaks makes it difficult 
to reliably determine the bubble's position. 

We present the results for R b = 20 Mpc in Figures 14.101 and 14.111 for the GMRT 
and MWA respectively. The HI fluctuations do not pose a problem for determining the 
position of such bubbles using the GMRT. In all the realizations of the GMRT simula- 
tions there is a peak at the expected position. The FWHM ~ 40Mpc is approximately 
the same along 9 x ,9 y and Uf and is comparable to the separation at which the overlap 
between the bubble and the filter falls to half the maximum value. The HI fluctuations 
does introduce spurious peaks, but these are quite separated from the actual peak and 
have a smaller height. We do not expect these to be of concern for position estimation. 

The MWA simulations all show a peak at the expected bubble position. The FWHM 
along 9 (~ 60 Mpc ) is somewhat broader than that along v (~ 40 Mpc). The low 
spatial resolution ~ 8 Mpc possibly contributes to increase the FWHM along 9. The 
HI fluctuations introduce spurious peaks whose heights are ~ 50 % of the height of the 
actual peak. 

4.3.4 Bubble detection in patchy reionization 

The SB scenario considered till now is the most optimistic scenario in which the HI 
traces the dark matter. The presence of ionized patches other than the one that we are 
trying to detect is expected to increase the contribution from HI fluctuations. We first 
consider the PR1 scenario where there are several additional ionized bubbles of radius 
6 Mpc in the FoV. Figures 14.121 & 14.131 show the mean value of the estimator (E) and 
3 — a error-bars as a function of Rf for the GMRT and the MWA respectively. These 



96 



4.3 Results 



-80 -40 40 80 -80 -40 40 80 




v, MHz v, MHz 



Figure 4.10: Same as the Figure KB for R b = 20Mpc for the GMRT. 
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Figure 4.12: The mean (E) and 3 — a error-bars of the estimator as a function of Rf for 
the GMRT estimated from 24 different realizations of the PR1 scenario. In 
all cases the filter is exactly matched to the bubble. 

were estimated from 24 different realizations of the simulation, using a filter exactly 
matched to the bubble. 

We find that the results are very similar to those for the SB scenario except that the 
signal is down by 0.6 due to the lower neutral fraction (xhi = 0.62) in the PR scenarios. 
Ionized bubbles with radius Rb = 8 Mpc and = 12 Mpc or smaller cannot be detected by 
the GMRT and MWA respectively due to the HI fluctuations. These limits are similar 
to those obtained in simulations of the SB scenario. 

In the PR2 scenario the FoV contains other ionized bubbles of the same size as the 
bubble that we are trying to detect. We find that bubble detection is not possible in 
such a situation, the HI fluctuations always overwhelm the signal. This result obviously 
depends on number of other bubbles in the FoV, and this is decided by xui which we 
take to be 0.62. A detection may be possible at higher xm where there would be fewer 
bubbles in the FoV. 
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Figure 4.13: Same as the Figure H?H for the MWA. 
In the SM scenario, the very large computation time restricts us from generating 
several realizations with central ionized regions of different sizes. Hence we are unable 
to study the restriction on bubble detection. We have only three realizations all of which 
have the same ionized region located at the center of the box. Based on these we find 
that the mean estimator (E) is ~ 30 times larger that the standard deviation due to HI 
fluctuations. The detection of a bubble of the size present in our simulation (Figure l4~TT) 
is not restricted by the HI fluctuations. We present size determination results in Figure 
H~T4"1 We see that the SNR peaks at R f = 42 Mpc and not at R f = 27 Mpc. Recall that 
in the 21-cm map (Figure |4~T|) we have visually identified the former as the bubble's outer 
radius which includes several small patchy ionized regions towards the periphery and the 
latter is the inner radius which encloses the completely ionized region. We see that the 
matched filter identifies the bubble's outer radius. To study the effect of non-sphericity 
we compare our results in Figure 14.141 with the predictions for a spherical bubble of 
radius Rb = 42 Mpc embedded in uniform HI with the same neutral fraction xhi = 0.5. 
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Figure 4.14: Same as Figure fl~6l for the SM scenario for the GMRT. The dotted, dashed 
dotted and dashed lines show results for three different realizations. To 
show the effect of non-sphericity, we compare these results with predictions 
for a spherical bubble of sized Rb = 42 Mpc embedded in uniform HI with 
neutral fraction xm = 0.5(solid line). The vertical line at Rf = 28 Mpc 
shows the radius up to which the bubble is fully ionized and the SNR 
follows the spherical predictions. The SNR peaks at Rf = 42 Mpc marked 
by another vertical line. 
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We find that our results for the SM scenario follow the spherical bubble prediction up 
to a filter size Rf = 28Mpc (marked with a vertical line in Figure H.14I) . corresponding 
to the bubble's inner radius which encloses a perfectly spherical ionized region. Beyond 
this, and upto the outer radius of 42 Mpc, the HI is not fully ionized. There are neutral 
patches which introduce deviations from spherical symmetry and cause the SNR to fall 
below the predictions of a spherical bubble beyond Rf = 28 Mpc. The deviations from 
sphericity also broadens the peak in the SNR relative to the predictions for a spherical 
bubble. 

Our results based on the SM scenario show that the matched filter technique works 
well for bubble detection and for determining the bubble's size even when there are devi- 
ations from sphericity. We obtain good estimates for the extents of both, the completely 
ionized region and the partially ionized region. For the SM scenario, Figure 14.151 shows 
how well the bubble's position can be determined in a blind search. We have followed 
the same method as described for the SB scenario in subsection 3.3. We see that the 
SNR peaks at the expected position. Further, as the bubble size is quite large > 27 Mpc 
there are no spurious peaks. 

4.4 Redshift Dependence 

Results shown so far are all at only one redshift z = 6. It would be interesting and useful 
to have predictions for higher redshifts. However, addressing this issue through direct 
computations at different redshifts would require considerable computation beyond the 
scope of this work. Since we find that the analytic predictions of the Chapter 3 are in 
good agreement with the simulations of the SB scenario, we use the analytic formalism 
to predict how different quantities are expected to scale with increasing z. 

The redshift dependence of some of the quantities like the system noise, the back- 
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Figure 4.15: Same as Figure |4~§1 for the SM scenario for the GMRT. The x-axis shows 
the comoving distance of the filter position from the center of the box. The 
three curves respectively show results for a search along three 9 X , 6 y and v 
axes. 
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ground 21-cm brightness I u , and the angular and frequency extent of a bubble of fixed 
comoving radius causes the SNR to decrease with increasing z. On the other hand the 
z dependence of the neutral fraction, the baseline distribution function and the effective 
antenna collecting area acts to increase the SNR at higher redshifts. We find that with 
increasing z both (E) and \J ((A£') 2 )hf decrease by nearly the same factor so that the 
restriction on bubble detection does not change significantly at higher redshift in the 
SB scenario. Assuming that the neutral fraction does not change with z, the SNR for 
bubble detection decreases with increasing redshift, the change depending on the bubble 
size. For example, for the GMRT at z — 10 the SNR decreases by a factor ~ 7 and 6 for 
bubbles of size Rb = 10 and 20 Mpc respectively. For the MWA this factor is 3 for both 
these bubble sizes. We expect a similar change in the SNR for the patchy reionization 
scenarios. The drop in SNR is slower for the MWA relative to the GMRT because the 
effective antenna collecting area of the MWA increases at higher redshifts. 

The SNR is directly proportional to the global neutral fraction xm which increases 
with z. The details of how xm and the HI fluctuations change with redshift depends on 
how reionization proceeds with time, an issue beyond the scope of this thesis. 



4.5 Summary 



We ha ve used a visibility-based formalism, introduced in 



Datta. Bharadwaj fc Choudhury 



(120071 ). to simulate the detection of spherical HII bubbles in redshifted 21 cm maps 
through a matched-filtering technique. The main aim of this work is to use simulations 
to quantifying the limitations for bubble detection arising from the HI fluctuations out- 
side the bubble. We have computed the results for two instruments, namely, the GMRT 
and the upcoming MWA. Our main conclusions are as follows: 
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In the case where the HI fluctuations outside the bubble trace the dark matter 
distribution (SB scenario), we find that bubbles with radius Rb = 6Mpc and 
= 12 Mpc or smaller cannot be detected by the GMRT and MWA respectively due 
to the HI fluctuations. Note that this limitation is fundamental to the observations 
and cannot be improved upon by increasing integration time. 

For targeted observations of ionized bubbles, the bubble size can be determined to 
an accuracy limited by the instrument's resolution; we find that HI fluctuations 
do not play any significant role. However, the HI fluctuations can restrict the 
accuracy to which the neutral fraction can be estimated. In addition, we find that 
determining the position of the bubble in a blind search could be quite difficult 
for small (~ 10 Mpc) bubbles as the HI fluctuations introduce large fluctuations 
in the signal; for larger bubbles the accuracy is determined by the instrument's 
resolution. 

In a scenario of patchy reionization where the targeted HII region is surrounded 
by many small ionized regions of size ~ 6 Mpc (PR1 scenario), the lower limit for 
bubble detection is similar to that in the SB scenario. Thus the assumption that 
the HI traces the dark matter gives a reasonable estimate of the contribution from 
HI fluctuations if the background ionized bubbles are small ~ 6 Mpc. However, the 
situation is quite different when the surrounding bubbles as of similar size as the 
targeted bubble (PR2 scenario). The large HI fluctuations do not permit bubble 
detection for a neutral fraction xhi < 0.6. Thus for xm = 0.6 or lower, bubble 
detection is possible only if the other ionized regions in the FoV are much smaller 
than the bubble that we are trying to detect. 

The matched filter technique works w ell for more realistic cases based on th e semi- 



numeric modelling of ionized regions (jChoudhury. Haehnelt fc Regan 



. Here 
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the bubbles are substantially non-spherical because of surrounding bubbles and 
inhomogeneous recombination. Our method gives a good estimate of the size of 
both the fully ionized and the partially ionized regions in the bubble. 



To put our conclusions in an overall perspective, let us consider an ionized bubble 
around a luminous OS O at z > 6. We expect Rh > 30Mpc from studies of QSO 



absorption spectra (I Wyithe &: Loeb 



2004 



Mesinger fc Haiman 



20041 ). It has also been 



pointed out that these bu bbles may survive as la rge "gray fossils" a long time after 



the source has shut down ( IFurlanetto et al 



20081 1 . It will be possible to detect such 



bubbles o nly if the background b ubble s are smaller, say, < 30 Mpc. We find from 



models of iMesinger &: Furlanettol (120071 ) that the typical sizes of ionized regions when 
xhi ~ 0.3(0.1) is ~ 20(70)Mpc. Though these values could be highly model-dependent, 
it still gives us an idea that the bubbles around the luminous QSO would be detectable 
even in a highly ionized IGM with, say, xm ~ 0.3. If the size of the targeted bubble is 
larger, then this constraint is less severe. This gives a realistic hope of detecting these 
bubbles at z > 6 with near-future facilities. 

A caveat underlying a large part of our analysis is the assumption that the bubbles 
under consideration are perfectly spherical. This is note the case in reality. For exam- 
ple, non- isotropic emission from the sources (QSOs), density fluctuations in the IGM 
and radiative transfer effects would distort the shape of the bubble. The semi-numeric 
simulations (SM scenario) incorporate some of these effects and give an estimate of 
the impact of the deviations from sphericity on bubble detection. This is an impor- 
tant issue which we plan to address in more detail in future. In addition, the finite 
light trav el time gives rise to an apparent non-sphericity even if the physical shape is 



spherical (IWyithe &: Loeb 



2004a : 



Yu 



20051 ). This effect can, in principle, be estimated 



analytically and incorporated in the filter. We plan to address this effect in future. 
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5 Optimum Redshift for Detecting 
Ionized BubblesL 



5.1 Introduction 

In Chapter 3 we have introduced a matched filter technique to detect individual ionized 
bubbles in 21-cm maps. The technique optimally combines the entire signal from an 
ionized bubble while minimizing the noise and foreground. In Chapter 4 we use different 
simulated 21-cm maps to study the impact of the HI fluctuations outside the bubble 
that we are trying to detect. 

Various redshift dependent parameters are important in bubble detection. The back- 
ground 21-cm brightness and the angular and frequency extent of a bubble of a fixed 
comoving size increase from higher to lower redshift. The system noise in low frequency 
radio experiments is dominated by the sky temperature which becomes lower at higher 

* A sharpened, compact version of this chapter is presented in the paper 
"The optimal redshift for de t ecting ionized bubbles in HI 21-cm maps" by 
\Datta. Bharadwai & Ch oudhurv. (2009 ) 
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frequency. This makes the noise lower at lower redshifts. All the above mentioned fac- 
tors are expected to enhance the signal to noise ratio (SNR) at lower redshifts. On the 
other hand, the neutral hydrogen fraction decreases as the reionization process proceed 
and it becomes almost zero around redshift z — 6. We will also see later in the text 
that the baseline distribution also changes significantly with the frequency band of our 
interest. This change is also expected to enhance the SNR at higher redshifts. Again, for 
MWA like experiments the increase in effective collecting area of the individual antenna 
is expected to reduce the system noise at higher redshifts. Together these parameters 
will try to suppress the SNR at lower redshifts with no signal at redshift z = 6 where 
the Universe is completely reionized. Because of these two oppositely behaving set of 
parameters there is an intermediate redshift where the SNR for a bubble of fixed co- 
moving size peaks. The detection of ionized bubbles is optimal at this redshift. The 
prior knowledge of this optimum redshift is important for any future attempt to detect 
ionized bubbles in 21-cm maps. 

In this Chapter we investigate the optimum redshift to detect ionized bubbles con- 
sidering different models of the neutral fraction evolution with redshift. Predictions are 
made for the GMRT and the MWA. We also establish scaling relations for the SNR with 
redshift for a fixed value of bubble size. 

The Chapter is organized as follows. We begin with a brief description on the 21-cm 
signal from an ionized bubbles and its redshift dependence Section 5.2. In Section 5.3 
we present a short note on our matched filtering technique to detect ionized bubbles. 
In Section 5.4 we establish scaling relations for the SNR with redshift considering both 
uniform-frequency independent and nonuniform-frequency dependent baseline distribu- 
tion. In the same section we also briefly describe the models of HI evolution that we 
adopt for our calculations. We present our results in Section 5.5 and conclusions in 
Section 5.6 
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5.2 Ionized Bubbles m Redshifted 21-cm Observations 

Throughout out this Chapter we adopt cosmological parameters from Dunkley et 
al.(2008). For the GMRT we use the antenna specifications from their website and for 
the MWA we use the instrumental parameters from Bowman et. al. (2007). 

5.2 Ionized Bubbles in Redshifted 21-cm Observations 

— * 

The quantity measured in radio-interferometric observations is the visibility V(U, v) 
which is measured in a number of frequency channels v across a frequency bandwidth 
B for every pair of antennas in the array. The visibility is related to the sky specific 

— * 

intensity I v {9) as 

v(u, u) = J d 2 eA(e)i v (ey^ ld (5.1) 

— * 

where the baseline U = d/A, d is physical separation between a pair of antennas pro- 
jected on the plane perpendicular to the line of sight. A and 9 respectively are observing 
wavelength and two dimensional vector in the plane of the sky with origin at the cen- 

— * 

ter of the field of view, and A{9) is the beam pattern of the individual antenna. We 
have used A(9) described in Chapter 3. The visibility recorded in radio-interferometric 
observations is a combination of four separate contributions 

V(U, v) = S(U, u) + HF(U, u) + N(U, u) + F(U, u) (5.2) 

where S(U, v) is the HI signal that we are interested in, HF(U, v) is contribution from 
fluctuating HI outside the target bubble, N(U, v) is the system noise which is inherent 
to the measurement and F(U, v) is the contribution from other astrophysical sources 
referred to as the foregrounds. 

We consider a spherical ionized bubble of comoving radius Rb centered at redshift 
z c at the center of the FoV. This bubble is also assumed to be embedded in an uni- 
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form inter-galactic medium (IGM) with neutral hydrogen fraction a; hi- The planar 
section through the bubble at a comoving distance r v is a disk of comoving radius 
R v = Rb\/l — {Av/ Av b ) 2 where Av = v c — v is the distance from the the bubble center 
v c in frequency space with v c = 1420MHz/(l + z c ) and Au b = Rb/r' Vc is the bubble size 
in the frequency space. The expected observed visibility from this bubble can be written 
as (for details see Chapter 3) 



^center (U , v) = -7lI u X m 9l 



2J 1 (2nU9 h 



2%U9 U 

where I v = 2.5 x 10 2 ^ ) (x) (^i)) is tlie rad iation background from the uniform 
HI distribution. 6 U = R u /r u is angular radius of the circular disc through the bubble 
at comoving distance r v . J\{x) and Q(x) are the first order Bessel function and the 
Heaviside step function respectively. 

The peak value of the signal is 5(0, v) = nl u 6l. For a fixed value of bubble size the 
peak value changes with redshift. The Bessel function J±(x) has the first zero crossing 
at x — 3.83. As a result, in baseline the signal S(U, v) extends to Uq = O.Gl^ 1 where it 
has the first zero crossing. The signal extends over Av = ±Aub in frequency. Note that 
the angular size 9 V and extent in frequency of a bubbles for a fixed value of comoving 
size change with redshift because of change in comoving distance r v and r' u . We will 
discuss the effect of these change on bubble detection at various redshifts in Section 5.4. 
Detecting these ionized bubbles will be a big challenge because the signal is buried in 
noise and foregrounds which are both considerably larger in amplitude. Whether we are 
able to detect the ionized bubbles or not depends critically on our ability to construct 
optimal filters which discriminate the signal from other contributions. 
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5.3 Matched Filtering Technique to Detect Ionized 
Bubbles 



For optimal signal detection we consider the estimator (Chapter 3) 



E 



Vb) 



a,b 



. o,,b 



(5.4) 



where Sf(U,u) is a filter which has been constructed to detect a particular ionized 
bubble, V(U a , u b ) respectively refer to the observed visibilities and U a and v h refer to 
the different baselines and frequency channels in the observation and in eq. (15.4ft we are 
to sum over all independent data points (visibilities). The filter Sf(U,is) depends on 
[Rf, z c , 6 C ] the comoving radius, redshift and angular position of the bubble that we are 
trying to detect. 

We now calculate (E) the expectation value of the estimator. Here the angular 
brackets denote an average with respect different realizations of the HI fluctuations, 
noise and foregrounds, all of which have been assumed to be random variables with zero 
mean. This gives (V(U, v)) = S(U, v) in the continuum limit and 



(E) 



d 2 U dvp N (U,v) S f *(U,v)S(U,v) 



(5.5) 



The variance of the estimator which is the sum of the contributions from the noise 
(NS), the foregrounds(FG) and the HI fluctuations (HF) can be written as 



((A^) = ((A^) ns + ((A^) fg+ ((A^) hf 



(5.6) 
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Assuming that the noise in different baselines and frequency channels are uncorrelated 
we have 



((AE)% S = a 2 j d 2 U jdup N (U,u) | S f (U,u) 

where a is the rms. noise in the image. 

The variance due to HI fluctuations is given as 



(5.7) 



(AE) 



2 ) 

/ HF 



xp N (U, vx)pn0, v 2 )S f *(U, vi)Sf{U, u 2 ) 

xC 27r u{Vl,V2) 



(5. 



where CW/(Vi, V'l) is the multi-frequenc y angu lar power spectrum of the HI fluctua- 



tion from 



Datta. Choudhury fc Bharadwaj 



(120071 ) (Chapter 2). 



5.4 Redshift Dependence: Scaling Relations 

5.4.1 Uniform and frequency independent baseline distribution 

The estimator (E) and its variance depend on the size (Rb), redshift (z c ) and position 
(9 C ) of the bubble. In this Section we present a scaling relations for the estimator 
describing its change with redshift. We assume a bubble of fixed size at the center of 
the FoV and the filter is exactly matched with bubble. In a situation where the baseline 
distribution (p(U,u)) is (i) uniform in baseline over the signal S(U,u) (ii) independent 
of frequency channels and (iii) the frequency bandwidth is much larger than the bubble 
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(BW » Au b ) the estimator scales exactly as 

(E)(z)«x 2 m (z)PXAv b (5.9) 

where xui(z) quantifies the evolution of the neutral hydrogen fraction. The redshifted 
21-cm background I v scales as (l+z)^ 1 - 5 at higher redshifts. The angular size 9 Uc which is 
oc \jr v approximately scales as (1 + z)~~ ' 25 for a fixed value of bubble size. Similarly the 
frequency extent of the bubble (2 Au b ) decreases with redshift as (l + z)~ - 5 . Considering 
all these quantities we get 

{E){z)kxUz){1 + z)-\ (5.10) 

The system noise is dominated by the sky temperature T s k y for the frequency range 
of interest. The sky temperature T sky increases with observed wavelength and scales 
as (1 + z) 13 , where (3 is the spectral index. Considering this fact and under the same 
assumptions mentioned above the noise variance can be written as 

((AE) 2 ) m (z)^x 2 m (z)(l + z)^-\ (5.11) 

In this equation we have assumed that the individual antenna collecting area is inde- 
pendent of frequency (applicable to the GMRT). For MWA like antenna the individual 
antenna collecting area is expected to increase with observed wavelength as ~ A 2 and 
the noise variance scales as 

((A£) 2 > NS (^)oc(l+z) 2 ' 3 - 8 . (5.12) 

The fluctuating HI outside the ionized bubble also contribute to the variance of the 
estimator. This contribution can not be reduced by increasing observation time and put 
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fundamental restrictions in detecting ionized bubbles. The HI distribution outside the 
target bubble during reionizaintion is highly unknown and differs drastically for different 
models of reionizations. For simplicity we assume that the HI outside the bubble trace 
the dark matter distribution. The presence of other ionized regions outside the bubble 
are expected to increase this contribution. The assumption that the HI outside the 
bubble trace the dark matter distribution will give the minimum variance contribution 
to the estimator. Analytical derivation of the scaling relation for the ((A.E) 2 )hf is not 
straight forward because of the term G%^\j{v\^ v<2). Assuming that the HI to trace the 
dark matter we numerically calculate ((A.E) 2 )hf at various redshifts using equation 15.81 
We then fit the data with a power law of (1 + z) which comes out to be 

((Ai) 2 ) HF oc4 I (z)(l + ^)- ia5 . (5.13) 

Detection of bubble is possible when (E)(z) > ((AE) 2 )np. When this is satisfied we are 
interested only in (E)/ J ((AE) 2 )ns which we call SNR. Note that the term xui(z) in the 
equations 15.91 15.111 is almost constant along the line of sight of the ionized bubble. This 
factor does not affect the value of the quantity SNR. One can drop x hi {z) and similarly 
I u from the filter Sf(U, v) without losing the effectiveness of the method. We calculate 
the SNR at different redshifts ranging from z = 6 to 16 for a fixed value of bubble size. 
Using equations 15.101 and 15.111 for a GMRT like antenna we have 



SNRoax m (z)(l + z)' 13 - 2 (5.14) 
For a MWA like antenna this is 

SNRocx m (z)(l + z)- 13 (5.15) 
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Note that in the above two equations xm(z) is a increasing function of redshift z and 
for a typical value of f3 m 2.5 — 2.8, we expect the SNR to be peaked at an intermediate 
redshift during reionization. As we discussed earlier this is due to three reasons, firstly 
the sky temperature increases at low frequency ie, at higher redshifts. Second the back- 
ground specific intensity l v also decreases at higher redshift. Finally the angular size 
and frequency extent of the bubble also become smaller at higher redshift due to the 
increase in r v and r' v with redshift . We see that the two factors are opposite in nature- 
one increases with redshift whereas the other decrease. There is an intermediate redshift 
z where the SNR is maximum. 

5.4.2 Non uniform and frequency dependent baseline distribution 

The baseline distribution in general is not uniform in U (see figure 3.5 in Chapter 3) 
and also frequency dependent. This non uniformity changes the above scaling relations. 
We introduce a parameter n to take into account this effect and now modified SNR can 
be written as 

SNR oc x m (z)(l + Z ) n ' 2 -P- 2 (5.16) 
for the GMRT like antenna. For MWA like antenna this is 

SNRocx m (z)(l + z) n/2 ~P (5.17) 

The normalized baseline distribution function Pn{U, v) is frequency dependent and the 
parameter n takes care of that. Note that n = is situation for uniform baseline 
coverage. The value of n depends on the bubble size. For large bubbles the signal 
remains confined within small baselines whereas for smaller bubbles a significant amount 
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of signal spreads to larger baselines. Thus the normalized baseline distribution function 
Pn{U, v) acts in different manner on bubbles of different sizes. For a situation when the 
Pn{U, v) is constant in U over the signal we have n = 2. In a situation where the signal 
is constant in U over the whole baseline range Pn{U, v) does not play significant role 
and n — 0. In general we have < n < 2. We will discuss more about the exact values 
of n for the GMRT and the MWA later. 

The GMRT baseline distribution is well described by the sum of a Gaussian (at small 
baselines) and an exponential (at large baseline) (Figure 3.5 in Chapter 3). Here Pn{U, v) 
remains almost constant up to a baseline U = 100 for the frequency range of interest and 
then it decays. For the larger bubbles of radius Rb > 50Mpc, the signal is confined with 
a maximum baseline of U — 100 and we have n = 2 for which SNR oc xm(z)(l + z)~ 3S 
where we assume (3 = 2.6. For smaller bubbles the signal extends to larger baselines 
where the baseline distribution function Pn{U, v) decays with U and the n value also 
becomes smaller. For example, for the bubbles of sizes Rb = 30, 20 and lOMpc the 
values of n are 1.27, 1.07 and 0.5 respectively. We have assumed that the MWA antenna 
distribution is p ant ~ r -2 , where r is physical distance from the center of the array. Here 
Pn(U, v) also decays monotonically. For the bubble sizes Rb = 50, 30, 20 and 10 Mpc the 
values of n are 0.45, 0.18 and respectively. 

The behavior of the HI fluctuations contribution also changes with redshift when 
we calculate for non-uniform and frequency dependent baseline distribution. As we go 
to higher redshifts the HI fluctuations decreases but slowly than uniform baseline case. 
This is due to the same reason as discussed in the 1st paragraph of this subsection. Using 
equation [531 we calculate ((AE) 2 )^ numerically at different redshifts for the bubble sizes 
R b = 10,20,30 and 50 Mpc. For the GMRT we find that ^ ((AE) 2 ) UF ~ x m (z)(l + 
z)~ 3 - 67 , xm(z)(l + z)~ 3 - 37 , xm(z)(l + z)~ 3 - 25 and xm{z){l + z)~ 3 - 16 respectively. For the 
MWA these are x m (z)(l + z)-^ 26 , x m (z)(l + z)~ AAS , x m (z)(l + z)- A and ^hi^XI + ^T 3 - 85 
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Figure 5.1: This shows the evolution of the mean neutral fraction xhi with redshift for 
the two different reionization models discussed in the text. 

respectively. 



5.4.3 Evolution of neutral fraction with redshift 



In this work, we consider two physically motivated models of reionization, namely, the 
early reionization model and the late reionization model which we call Model I and 



Model II respectively. T 



rtChoudhury fc Ferrara 



rese models are constructed using t 



2005 



Choudhury fc Ferrara 



re semi-analytical formalism 



2006al ) which implements most of 



the relevant physics governing the thermal and ionization history of the IGM, such as 
the inhomogeneous IGM density distribution, three different classes of ionizing photon 
sources (massive PopIII stars, PopII stars and QSOs), radiative feedback inhibiting 
star formation in low-mass galaxies and chemical feedback for transition from PopIII 
to PopII stars. The models are consistent with various observational data, namely, the 
redshift evolution of Lyman-limit absorption systems, the Gunn- Peterson effect, electron 
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scattering optical depths, temperature of the IGM and cosmic star formation history. 

In Model I, hydrogen reionization starts around z ~ 16 driven by metal-free (PopIII) 
stars, and it is 50% complete by z ~ 10. The contribution of PopIII stars decrease below 
this redshift because of the combined action of radiative and chemical feedback. As a 
result, reionization is extended considerably completing only at z ~ 6 (Figure I5TT1) . 

In Model II, the contribution from the metal-free stars is ignored, which makes reion- 
ization start much later and is only 50 per cent complete only around z m 7.5. The main 
difference between this model with the previous one is in their predictions for the electron 
scattering optical depth (which is 0.15 for the early reionization model and 0.06 for the 
late reionization m odel). Recent measu rements of the electron scattering optical depth 



r e = 0.087 ± 0.017 flDunklev et al. 



20081 ) suggests that these two models are possibly at 



two extreme ends of different possible reionization scenarios. 



5.5 Optimum Redshift to Detect Ionized Bubbles 

We consider bubbles of comoving size ranging from 2 Mpc to 50 Mpc at redshifts from 
z = 6 to 16. Figure 15.21 shows the SNR contours for 1000 hrs of observations with 
the GMRT. The left panel shows results for a constant neutral fraction xm = 1- This 
shows the joint effect of factors such as the background 21-cm brightness I u , angular 
and frequency extent of the bubble, sky temperature, normalized baseline distribution 
function and the effective area of individual antenna. We see that the SNR is maximum 
at redshift the z = 6 for any bubble size. As we have discussed in Section 5.4 (equations 
I5.16l and l5.17p this is because of the lower noise at lower redshift. The region to the left 
of a line is allowed for the detection. For example a 5 a detection is possible for a bubble 
size Rb > 16 Mpc at redshift z = 6. At z > 11 a 5cr is not possible even for bubbles of 
size 50 Mpc. The middle shows the same results for Model I. We see that the redshift 
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Figure 5.2: This shows the SNR contours for 1000 hrs observations with the GMRT. The 
left panel shows results for a constant neutral fraction xm = 1. Middle and 
right panel show results for the Model I and the Model II respectively. The 
HI fluctuations dominate over the signal in the shaded region and bubble 
detection is not possible. 




Figure 5.3: Same as the figure 15721 for the MWA. 
range z ~ 9 to 10 where the SNR is maximum is the optimum redshift range for bubble 
detection. We also see that a 3 a detection is possible for bubbles of size > 46 Mpc at 
the optimum redshift. A 3 a detection is also possible for a bubble size Rb > 50 Mpc 
in the redshift range z ~ 7 to 11 . Results for the Model II are presented in the right 
panel. Here we see that redshift z = 8 is optimum for bubble detection. In this model 
the reionization occurred at lower redshift and hence the SNR is comparatively higher. 
Bubble detection is relatively easier in Model II. For example a 5a detection is possible 
for Rb > 38 Mpc at z = 8 which is not possible in the Model I. In all the panels bubble 
detection is not possible on the shaded region due to the HI fluctuations. 

Figure 15.31 shows the same results for the MWA. The effective collecting area of the 
individual MWA antenna is expected to increase with wavelength i.e, with redshift as 
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Figure 5.4: This shows the SNR contours for a fixed value of bubble size with the GMRT. 

Left panel shows results for a Rb = 20 Mpc bubble and 4000 hrs of observa- 
tions. Right panel shows results for Rf, = 50 Mpc bubble size and 1000 hrs 
of observations 

~ (1 + z) 2 . This will reduce the noise rms at higher redshifts and make bubble detection 
easier at higher redshifts than the GMRT. For example a 5a detection is possible in 
Model I for bubble Rb = 50 Mpc in the redshift range z = 9 to 15 which seems impossible 
with the GMRT. The SNR is 3 times higher than the GMRT at redshift z = 16 for a 
Rb = 46 Mpc bubble. The optimum redshift in Model I for the MWA is z ~ 11 which is 
slightly higher than the GMRT. Detection of bubbles in model II is more optimistic in 
the MWA where 5a is possible for bubbles of Rb > 38 Mpc at redshift z ~ 8.5 which is 
the optimum redshift of detection in this model. 

In the above discussion we consider different models of the mean HI evolution. The 
HI evolution surrounding the ionized bubbles is highly unknown and model dependent. 
Keeping this in mind we consider a bubble of a fixed size and calculate the SNR for 
neutral fraction xhi = to 1 at various redshift ranging from z = 6 to 16. We present 
results in Figures 15.41 and 15.51 for the GMRT and the MWA respectively for bubble size 
Rb = 20 Mpc (left panel) and 50 Mpc (right panel). The left panel in the figure I5~4l 
shows results for 4000 hrs of observation. The Figure demonstrates the detectibility of 
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Figure 5.5: Same as the figure [53] for the MWA. 
a Rb = 20(50) Mpc bubble embedded in the HI of some neutral fraction. We see that a 
higher neutral fraction is required to detect bubbles at higher redshifts. For example we 
see that a 3<r detection for a bubble Rb = 50 Mpc is possible for neutral fraction ~ 0.1 
at redshift z = 6 with the GMRT (Figure 15.41) but the neutral fraction ~ 1 is required 
at z = 12 to detect the same bubble. Bubble detection is found to be easier at higher 
redshifts for the MWA. 

5.6 Conclusions 

We investigate the optimum redshift for bubble detection. We find that for early reion- 
ization the redshift z = 9 and 11 are respectively the optimum redshifts for the GMRT 
and the MWA. For late reionization redshift z = 8 is found to be the optimum redshift 
for both the experiments. Bubble detection will be easier if the reionization occurred 
late. This is because of less noise and higher 21-cm background intensity. The MWA is 
suitable for bubble detection at higher redshifts (z > 8.5) than the GMRT. We also find 
that at redshift z — 6 if surrounding mean neutral fraction is xm > 0.1 a bubble of size 
Rb > 50 Mpc is possible to detect with 1000 hrs of observation for both the experiments. 
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Appendix 



A Calculation of the angular power spectrum C\{y\,ui) 

In this section, we present the details of the calculation for the 21 cm angular power 
spectrum C\(y\, 1*2). The first step would be to calculate the spherical harmonic compo- 
nent ai m of the brightness temperature T(u, n). Using the expression (12. 6p for f/ HI (k) in 
equation (12.101) . the expression for a/ m can be written as 



f xm / dn Y*(h 



d 3 k 



x 



A m (k) + (n • k) 2 A(k) 



-ifcr„(k-n) 



(Al) 



which then essentially involves solving angular integrals of the forms J dn Y lm e _lfcrv ( k ' n ) 
and j dfl Y lm (n • k) 2 e _lA:rv ( k ' n ) respectively. Expanding the term e - lfc M k ' n ) i n terms of 
spherical Bessel functions ji(kr u ), one can show that 



dn Y lm (h) e - ifer "( k - fl ) 



M-i) l ji(kr u )Y lm (k) 



(A2) 



Differentiating the above equation with respect to kr v twice 



dn (k • n) 2 Y^(n) e 



-ifcr„(k-n) 



-47r(-i) ( j'{{kr v )Y{ m {\i) 



(A3) 



where j"(x) is the second derivative of ji(x) with respect to its argument, and can be 
obtained through the recursion relation 



(2/ + l)fi\x)J-^-^-jUx) 



I 2 ^ (/ + l) 2 



21-1 21 + 3 
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21 + 3 

So the final expression of ai m is given by 



(1 + l)(/ + 2) ... 
^ J) o 3i+2{x) (A4) 



/d 6 k 

x [A H i(k) Ji (A;r ! ,)-A(k)j; / (A:r i ,)] (A5) 

The next step is to calculate the the power spectrum Ci(h>i,v 2 ) = (a/ m (^i) a z* m ( I/ 2))- 
The corresponding expression is then given by 

Ci{v u u 2 ) = {4n) 2 f ( Zl )f (z 2 ) x m ( Zl )x m (z 2 ) 

x ([A m (z u kMhr^) - A(z u kOjf^rvJ] 
x [A^ I (z 2 ,k 2 )^(A; 2 r i , 2 ) - A*(z 2 , ^)ji{k 2 r V2 )}) 

(A6) 

where we have put back the redshift-dependence into the expressions for clarity. Now 
note that we would mostly be interested in cases where v 2 — v\ = Au <C v\. In such 
cases, one can safely assume T(z 2 ) m T(z 1 ) and x m (z 2 ) m xhi(^i). Furthermore, 
the terms involving the ensemble averages of the form (AA*) can be approximated 
as (A(zi, \s. 1 )A*(z 2 , k 2 )) ~ (2-K) 3 5 D (k 1 — k 2 )P(z 1 , k{) and similarly for terms involving 
A H i- We can then use the Dirac delta function 5 D {ki — k 2 ) to compute the k 2 -integral, 
and thus can write the angular power spectrum as 

C,(Ai/)=G(i/,i/ + Ai/) 
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B Correspondence between all-sky and Hat-sky power spectra 



(4vr) 2 T 2 4i/ 




Y; m (k)Y lm (k) 



x Ji(kr u )ji{kr U2 )P A 2 (k) 



{ji{kru)ji{ kr u 2 ) + 3i(kr„ 2 )j" (kr u )}P Am (k) 



+j' l '(kr v )j'/(kr U2 )P(k)} 



(A7) 



Using the normalization property of the spherical harmonics J dn|Yi m (n)| 2 = 1, one can 
carry out the angular integrals in the above expression, and hence obtain the final result 
(12.111) as quoted in the main text. 

B Correspondence between all-sky and flat-sky power 
spectra 

As discussed in section 12.2.31 we shall mostly be interested in very small angular scales, 
which corresponds to i 3> 1. For high values of /, it is most useful to work in the flat-sky 
approximation, where a small portion of the sky can be approximated by a plane. Then 
the unit vector n towards the direction of observation can be decomposed into n = m+9, 
where m is a vector towards the center of the field of view and 9 is a two-dimensional 
vector in the plane of the sky. 

Without loss of generality, let us now consider a small region around the pole 9 — > 
0. In that case the vector 9 can be treated as a Cartesian vector with components 
{9 cos 0, 9 sin <$}. This holds true for any two-dimensional vector on the sky, in particular 
U = {U cos (pu, U sin0[/}. Then the spherical harmonic components of T(V, n) [defined 
in equation (12.101) ] can be written as 




(Bl) 
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where we have replaced J dQ — > J d8. Now use the expansion 

e - 2 ^u.e = J2(-i) m J m {2TiUey m{<t>u -^ (B2) 

m 

where J m (x) is the ordinary Bessel function. Further, we use the approximation for 
spherical harmonics 

Yi m (9, <P) J m (l0)^^ (B3) 

to write 

e -™™ « ^1 X)(-i) TO ^,m(^ ^) e_im ^ ( B4 ) 

m 

Then the two-dimensional Fourier transform of the brightness temperature [defined in 
equation (I2.16p ] will be 



f(i/,U)= y dfl e - 2 " iu e T(z/,n) 



m -1IU0U 



^ EH) 

_L ^ ( _i)m e -im^ ajhr ^ m(l/ ) (B5 ) 



m 

where we have used the expression (IB II) for a^ m in the last part. This gives a relation 
between the flat-sky Fourier transform T(z/, U) and and its the full-sky equivalent ai m {v). 
Using the above relation, we can calculate the power spectrum 



(TKU)T> 2 ,U0>~^/-^ ^(-i)"' e - im ^ e im '"" 



x(a2^, m (^i)aLc/',m'( z/ 2)) (B6) 
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B Correspondence between all-sky and Hat-sky power spectra 



Use the definition (a; m (z/ 1 )af, m ,(f 2 )) = Ci^u'^mm' and the property 

J2e--M<t>v-M = 2^\<Pu - M (B7) 

m 

to obtain 

(TKU)T> 2 ,U')> = 2tt C^ v {v u v 2 ) ^ 5 { S{4>u- M (B8) 

The last step involves writing the right hand side of the above equation in terms of the 
two-dimensional Dirac delta function, which follows from the expansion 

4 2) (U - U') = J d9 e-W- u 'H (B9) 

The exponentials can be written in terms of the spherical harmonics using equation 
(El: 

mm' 

(BIO) 

Finally use the orthonormality property of spherical harmonics J dOY^O, <ft)Yi' m >(6, <j)) = 
Sii'Smm' and the relation (IB7I) to obtain 

4 2) (U - U') = 2vr ^Snfa - M (Bll) 

Putting the above relation into flB8p . we obtain equation (12.201) used in the final text. 
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C Relation between visibility-visibility correlation and 
MAPS 



In this appendix we give the calculations for expressing the two visibility correlation 
in terms of the Multi- frequency angular power spectrum (MAPS). We can write the 
visibility V(U, v) as a two-dimensional Fourier transform of the brightness temperature 
T(9, v) [see equation (13. ip ] 

V(U, v) = (^j J d 2 9A(6, v)T(6, u)e^ (CI) 

where (dB/dT) u is the conversion factor from temperature to specific intensity and 
A(6, v) is the beam pattern of the individual antenna. The visibility- visibility correlation 
is then given by 

(vA-Ovft^© ( 8B 



x (T(9, v x )T$, V2 )) e ^M-u 2 ) 



(C2) 



The correlation function for the temperature fluctuations on the sky would simply be 
the two-dimensional Fourier transform of the MAPS C^u^x, 1*2) 

(T{6, vjTfi, u 2 )) = j d 2 U C 2nU (vi, P2)e- 2m ^- t} (C3) 
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C Relation between visibility-visibility correlation and MAPb 



Using the above equation in equation in (102 j) . we obtain 



27 



xA{U 1 -U,v 1 )A{U 2 + -U^ 2 ) (C4) 

where A(U, v) is the Fourier transform of the beam pattern A{9, v). If the beam pattern 
is assumed to be Gaussian A(9, the Fourier transform too is given by a 

Gaussian function 

A(U, u) = n9 2 e- n2u2e " (C5) 
Hence, the visibility correlation becomes 

x J d 2 U C 27rt /(ivi,/v 2 ) 

xe -n 2 [(U 1 -U) 2 9 2 1 + (U 2 +U) 2 9l] /QgX 

where 9\ and 9 2 are the values of 9q at ^ and v 2 respectively. Now, since the two 
Gaussian functions in the above equation is peaked around different values of U, the 
integrand will have a non-zero contribution only when \Ui + U 2 \ < (it max[0i, 9 2 \)~ l . 
In case the typical baselines are much larger than the quantity (tt max[9\, 9 2 ])~ 1 , the 
integral above can be well approximated as being non-zero only when JJ\ = —U 2 . Then 

(V{U x ,v x )V0 2 ,v 2 ))*5^ (|f)^ 
xit 2 9 2 1 9 2 1 C 27tUi {v 1 ,v 2 ) 



139 



5 Appendix 



_ x (_m_\(d_B\ (d_B_\ 

(C7) 

which is what has been used in equation (13. 9|) . 

In the continuum limit, the Gaussian A(U, v) can be approximated by a delta func- 
tion, i.e., A(U,u) ~ 5jj(U) (which corresponds to the limit 6> — > oo); the visibility- 
visibility correlation is then given as 

(V(U 1 ,v 1 )V0 2 ,v 2 )}=6%\u 1 + U 2 ) (^pj {^j 

xCtorUifaVi) (C8) 

which corresponds to equation (13.211) in the main text. 
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